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DEVELOPMENT OF AN ORTHOTROPIC HOLE ELEMENT 


by 


C. V. Smith* 

J. W. Markham 

J. W. Kelley 

K. Kathiresan 


SUMMARY 

The objective of this contract was to develop and implement a finite element which 
adequately represents the state of stress in the region around a circular hole in 
orthotropic material experiencing reasonably general loading. This has been achieved 
through a complementary virtual work formulation of the stiffness and stress matrices 
for a square element with center circular hole/ with implementation in the NASTRAN 
computer program. This report discusses the theoretical foundations of the develop- 
ment, describes the use of NASTRAN to obtain solutions, and presents sample 
problem results to demonstrate the performance of the hole element. 


Consultant, Associate Professor, School of Aerospace Engineering, 
Georgia Institute of Technology. This report carries no endorsement 
by Georgia Institute of Technology. 
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INTRODUCTION 


In present flight vehicle structures, and for some time into the future, mechanical 
fastening is expected to be the most used joining methed, so that connection details 
contain many circular holes. These connections are subjected to very general external 
loadings, with the additional complications of load transfer through fasteners in the holes 
and increasing use of nonisotropic material in the connection. The resulting local stress 
state in the neighborhood of each hole is very complex, with stress levels and gradients 
higher than those found in most regions away from the holes. Despite these complexities, 
the high stresses require that an accurate evaluation of the local stresses around each 
hole be obtained to achieve an adequate estimate of performance of the joint, as 
measured by fatigue and fracture mechanics analyses. As usual, the analyst is in the 
unenviable position of needing the best results in the regions of greatest complexity. 

One method to determine the stresses, through analytical solutions, would appear to 
give those desired results; but this procedure can be discarded for all practical purposes. 
The list of available solutions is not adequate, and there is little hope that the list can 
ever be expanded to include all possibilities of interest. 

The finite element method is a very effective way to treat problems with complex 
geometry, material, and loading, but the use of conventional finite elements is not 
efficient in applications involving fastener holes. Because of the complexity of 
the stress state and the high stress gradients, an excessive number of elements must be 
used in the neighborhood of each hole. This increases the cost of even one solution 
for a specified geometry and loading, and makes it very expensive to perform the 
multiple solutions with varying geometry and loading which are normally required in 
design. 

What is required is an efficient special purpose finite element which provides a 
sufficiently accurate representation of local stresses without the need for many 


X 



conventional elements. This has previously been done for isotropic material 
([8], [9], [10] ), and this report presents the extension to orthotropic material. 
Complementary virtual work concepts are used to develop the stiffness and stress 
matrices for a square element with a center circular hole. The report presents the 
necessary theoretical background for the development and then presents sample 
problem results obtained through implementation into the NASTRAN computer 
program . 


* Numbers in brackets denote reference numbers. 
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SECTION 1 

THEORETICAL BASIS FOR DERIVATION OF THE STIFFNESS MATRIX 


1 . 1 STATEMENT OP THE PROBLEM 

The finite element is considered to be a thin plate in a state of plane stress with 
displacements on the external boundary specified so as to insure displacement con- 
tinuity with adjacent finite elements (see Figure 1). Conditions on the inner boundary 
will be specified in a later section. 

The variables of primary interest are the stresses because these quantities are 
needed for studies in fatigue and fracture. Therefore, it was decided to seek a solu- 
tion to the problem through a complementary work formulation which gives the desired 
stresses directly with no need for differentiations of displacements with subsequent loss 
of accuracy. 

1 .2 STATEMENT OF‘THE PRINCIPLE OF COMPLEMENTARY VIRTUAL WORK 
The complementary virtual work (CVW) can be expressed as follows: 

CVW = t y* |^6rj|€|dA - ij |6fJ |u| dA- |6tJ |u| 


where t = thickness of hole element, assumed constant 
A = area of hole element 

= that portion of the boundary over which displacements are specified, 
consisting of the outer boundary plus portions of the inner boundary 
as appropriate 

= that portion of the boundary over which tractions are specified, 

consisting of portions of the inner boundary as appropriate 

icl = U e y I ^ = column of strains in the real state 
I j L X y xyj 

•|u| = |u^ *^yj^ ~ column of displacements in the real state 

|u|’ = specified displacements on the boundary 

|6rl= Ist 6t 6t J= row of stresses in the virtual state 
^ X y xy 
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Figure 1 . Finite Element 
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|aFj=|_w 5 fJ = 


row of body forces in the virtual state 



row of surface tractions in the virtual state, acting 
on element of boundary with normal n. 


The virtual state must satisfy the following admissibility requirements: 

a) stresses and body forces satisfy the equations of equilibrium in the 
finite element domain: 


96t 36t 

^+_^+6F =0 

3x ay X 

a6r a6r 

=0 

ox 0 y y 

6t - 6t =0 

xy yx 


(2) 


b) stresses and surface tractions satisfy the equations 


6T - 6t n +5 t n 
X XX yx y 

6T “ 6t n + 6t n 
y xy X y y 


(3) 


where n^, n^ = components of unit outer normal at a point on the plate boundary. 

The real state must meet the following requirements: 

a) the strains and stresses are related by constitutive equations which can be 
written in the form 

where |t\= |t T T I ^ = column of stresses in the real state 
I ) Lx y xyj 

|cj= array of material coefficients 


( 4 ) 
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b) 


the stresses must satisfy equations of equilibrium in the finite element 
domain: 


ar 

ar 

X 

o 

II 

+ 

hx 

ay 


ar 

xy 

+ y -0 

ax 

ay 

r 

- T =0 

xy 

yx 


with body forces neglected. 


( 5 ) 


c) on the boundary , the stresses and the specified tractions must be 
related as follows: 


T - T n + T n 
X XX yx y 

n 

T = T n + T n 

y xy X y y 


( 6 ) 


Note that all admissible real states can be described as stress states which satisfy 
equilibrium in the domain and traction boundary conditions. 


The principle of complementary virtual work can be stated as follows. If a particular 
real state can be found for which the complementary virtual work vanishes for every 
admissible virtual state, then that particular real state will also satisfy compatibility 
requirements in the interior and displacement boundary conditions. Therefore, this 
particular real state satisfies all the requirements for a solution to the given problem 
and is, indeed, the solution. 


The expression for complementary virtual work, as given in Equation (1), can be 
simplified by considering only those virtual states for which 


and 


Uf\ = LoJ 
LstJ = LoJ 


in the domain 
on the boundary S^. 


4 



These additional admissibility requirements are not necessary, but they do serve to 
reduce the complexity of the solution process. With these simplifications, the comple- 
mentary virtual work can be calculated as follows: 


CVW = 


J u 


which is the form to be used in the derivations which follow. 


( 7 ) 


For future use, note from Equation (7) that if real displacements |u| exist for which 


strains <€> must be identically zero in the domain. 


0 for every admissible 


L«tJ, 


then it follows that the associated real 


1 .3 USE OF COMPLEMENTARY VIRTUAL WORK IN APPROXIMATIONS 

For problems with no available exact solution, the principle of complementary virtual 
work can be used to achieve an approximate solution. For the real state, begin with 
a family of admissible stresses in the form 


|r(x, y)| = |Tp(x, y)| + [a (x, y)] (8) 

3x1 3x1 3xNNxl 


where )T \ denotes a stress state which satisfies the 

I P/ 

equilibrium equations (including body forces 
if applicable) and all nonhomogeneous boundary 
conditions on S 

T 


and 



denotes a column of N as-yet unknown coefficients 

denotes a family of stress states which satisfy 
equilibrium equations with no body forces and 
homogeneous boundary conditions on all of S^. ■ 
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Therefore, each column of [a] represents on admissible stress state; and Equation (8) 
represents an admissible stress state for arbitrary values of the N coefficients . 

Note that any known features of the real state stresses, such as satisfaction of symmetry 
conditions, can be included in the functions which constitute the elements of array 
[a (x, y)] . 


The real strains are now determined by substituting into Equation (4). 


(9) 


To complete the specification of the real state, the displacements on are assumed 
in the form 

|u (x, y)| = [C (x/ y)] |q| (10) 

2x1 2 X F Fxl 


where |q|denotes the generalized displacements on S^; there are F degrees of 

freedom. Therefore, the i*^^ column of C represents boundary displacements 

th ^ 

associated with a unit value of the i degree of freedom; and by assumption. 
Equation (10) is a sufficiently accurate representation of the specified |u|on 


The virtual state is specified by 


1 6r} = [a (x, y)] I a I 
3x1 3x N N X 1 


If Equation (3) is written in the form 


{6t}=[n] {6r} 

2x1 2x3 3x1 


then 


or 


(6t}= [n] [a] {□} 

{6t}= [7j {o} 


( 11 ) 


( 12 ) 
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Everything is now known for substitution into Equation (7) for the complementary 
virtual work. 



Define ^ [^] 

Nxl Nx3 3x3 


3 X 1 




Nx3 


[c] [aJcIA 

3x3 3x N 


[G] = t/ [r7 
Nx F N X 2 


‘c ] ds 
2 X F 


Therefore, CVW - 




(13) 


(14) 


(15) 


(16) 


Now, within the approximation of Equation (11), every virtual state is described 
by a set of values for coefficients |a| . Therefore, if the CVW is to vanish for 
every admissible virtual state, it follows that 

The values of ^ are determined from 
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With |iS| known from Equation (18), the solution for the real stress state is found 
from substitution into Equation (8). 

Note that when the- material constant matrix |cj is positive definite, then in the 
real state 

W H ^ every |t| . 

If the stress state is chosen as|r| = follows that 

L^J[a]t[c][a]{4^ 0 for every Ijsl . 

Therefore, the integrand in Equation (14) is positive definite, which means that the 
array |^hJ is positive definite. This guarantees the existence of |^hJ \ so that the 
solution in Equation (18) exists. 


1 .4 OBSERVATIONS ON THE ARRAY G 

In Equation (11), it was not necessary to assume that the virtual stress state |6 t| has 
the same spatial distribution |^A (x, y)j as was assumed for the homogeneous real stress 
state. Furthermore, it was not necessary to assume that the number of coefficients |a| 
equals the number of coefficients in the real state. However, if the CVW is expressed 
by the simplified form in Equation (7), then the number of coefficients lal must equal the 
number of coefficients |/3| . Also, the virtual state stresses must satisfy homogeneous 
equilibrium equations and give zero traction values on the boundary S^; but these are 
exactly the admissibility requirements for the homogeneous real stresses, so the virtual 
and homogeneous real stresses might as well be assumed to have the same spatial distribution. 


The external work, 6Wg, done by the virtual state surface tractions moving through the 
real state boundary displacements is given by 


6W3 = 



(19) 
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where there Is no contribution from because [6lJ -|_ojon . Substitute Equations (10) 
and (12) into Equation (19): 


6w^ = [aJ [r]’’[c]* {q} 


( 20 ) 


Then substitute Equation (15) to conclude 

6W, = [=J [g] 

From Equation (20) it follows that the G.. element in array G (the element in the 
i^” row and j ” column) represents the boundary work done by the virtual tractions 
associated with coefficient a. moving through the real displacements associated with 
q^. Array JgJ will be called the boundary work array. 

Because the virtual state stresses satisfy homogeneous equations of equilibrium, it 

follows that the virtual boundary tractions evaluated from Equation (3) will be self- 

equilibrating. Furthermore, because the virtual stresses satisfy homogeneous traction 

boundary conditions on S , it follows that the non-zero virtual tractions on S must 
'IT ^ 

be self-equilibrating. Therefore, the external work 6Wg , as expressed in Equation (20), 
should be zero if and only if the real displacements |q| represent a rigid body motion of 
the boundary; and this must be true for every admissible virtual stress state (for 
every LaJ )• These requirements provide checks on the formation of the |gJ array and 
provide guidance in the selection of the arrays I^tJ and which are used to derive [g] 
through Equation (15). 


1 .5 REQUIREMENTS ON THE ASSUMED BOUNDARY DISPLACEMENTS 


Equation (10) expresses the displacements |u| at any point on the boundary S^ in terms 
of the specified degrees of freedom |q| . If the external work 6Wg is to vanish for a 
rigid body motion, then obviously there is a requirement that Equation (10) must include 
rigid body motion. Let|q||^g denote rigid body generalized displacements of S^. 

Then substituted into Equation (10), the resulting |u|must represent 

rigid body motion everywhere on S^. That is. 
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{u(x,y)}„B-[c(x, yjl {q}|^B 

which places requirements on the elements in the array j^^J. 

1 .6 REQUIREMENTS ON THE NUMBER OF STRESS COEFFiC lENTS 

Because the work 6Wg/ as given in Equation (20), must vanish for every virtual stress 
state moving through real state rigid body motions, it follows that as a necessary 
condition 

[g] {i}={o} (22) 

NxF Fxl Nxl 

should be satisfied only by displacements which are derived from true rigid body 
motions of the boundary S^. That is. Equation (22) should be satisfied by rigid body 
|q| and should not be satisfied by any non-rigid body|q|-. 

Equation (15) shows that the array has N rows and F columns, where N is the 
number of stress coefficients and F is the number of boundary degrees of freedom 
in Equation (10). Let R denote the number of rigid body degrees of freedom which 
are consistent with any special features (such as symmetry conditions) included in the 
elements of array |^C(x, y)J and which cause non-zero displacements on S^. (Note 
that there might be a rigid body degree of freedom which causes zero displacements 
on S^. The resulting |q||^g =|oj-will certainly give 6Wg = 0 from Equation (20) but 
this provides very little information as a check on the array joj .) 

Now, because there exist only R true rigid body motions, it follows that the most 
general nontrivial solution for|q| in Equation (22) must contain R independent 
constants. If the R constants are selected to be R components of|q|, it follows that 
it must be possible to solve Equation (22) for the remaining (F-R) components in terms 
of the R rigid body components. From this it follows that the rank of jcj , denoted 
by rank joj , must be equal to (F-R). 
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If rank |gJ > (F-^z then it would be possible to solve Equation (22) for more than 
(F-^ components expressed in terms of less than R independent constants. This 
would mean fewer than R rigid body motions, which means that some of the true 
rigid body modes have been suppressed. If the structure were to experience the 
suppressed rigid body motion, there would be non-zero 6Wg and non-zero stresses 
from Equation (18). This possibility will not occur with self-equilibrating tractions. 


If rank |Gj < (F-R), then it would be possible to solve Equation (22) for less than (F-R) 
components expressed in terms of more than R independent constants. This would 
mean more than R rigid body motions, which means that spurious rigid body modes 
have been introduced into the system, if the structure experiences a spurious rigid 
body motion, there will be zero 6Wg and zero stresses even though the imposed dis- 
placements are actually not rigid body. This is bad and should be avoided. 

The requirement that rank |gJ = (F-R) imposes the obvious requirement that (F-R), 
for if N < (F-"R) then it is impossible to have rank [g] = (F-R) . Note that there is no 
theoretical objection to N > (F-7) provided that all N equations in Equation (22) 
are satisfied by the rigid body displacements (this is guaranteed since the virtual 

stress states are self-equilibrating). 


In order to develop additional requirements on the choice of stress coefficients, consider 
the array I q 1 in the partitioned form 


hRB 

1 X F lx (F-S) 1 X S 


(23) 


where q now denotes any set of generalized displacements on S selected so that 
. L RBJ ^ ^ 

q = |0 is sufficient to suppress all true rigid body motion. Clearly, there is a 

L '^^-1 - _ . I - I . 

requirement that S 5 R; but the choice of elements in the array is not unique. 

Consistent with the partitioning of [qj , partition array [g] in the form 


["]= [ 

N X F N X (F-S) 


'1 


^RB 
N X S 




(24) 
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Then W,=[_aJ[G,]{q,}+|_aj[G,J{iJ; 

and if the true rigid body motions ore suppressed by setting = |_^J' 

external work 6Wg should vanish. 


6W,=L->J[g,]{<I,} =0 

This must be true for every [oj , which requires 

[<^,] = W ' 

N X (F-S) (F-S) X 1 N X 1 


(25) 


Because true rigid body motion has been suppressed/ Equation (25) should be satisfied 
only by the trivial solution |q • Therefore/ it follows that the rank of array 

must equal to (F-S)/ for suppose rank [G^j< (F-S); say rank = M. Then 
it would be possible to consider a virtual stress state consisting of the proper choice 
of M independent stress states/ denoted by |°ij/ such that 


6W3 = 


L°lJ [ ^l] { '’l} 

1 X M M X (F-S) (F-S) X 1 


would vanish for arbitrdry |a J and nontrivial |q ^ j- . This indicates existence of a 
spurious rigid body motion/ which is to be avoided. 


In summary/ the number of independent stress states in the virtual stress state must 
first be sufficient to give N ^ (F-R). Then it is necessary to insure that rank joj =(F-^. 
Finally/ it is necessary to insure rankjo^ = (F-S)/ and this must be true for every choice 
ofj^qj^gJ in Equation (23). 
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1 .7 TREATMENT OF BOUNDARY CONDITIONS 


As described earlier/ the domain boundary consists of two parts - an portion over 
which displacements are specified and an portion over which tractions are specified. 
In a complementary virtual work formulation, the requirements on the boundary. 
Equation (6), impose requirements on the admissible stress states. Equation (8). If these 
requirements can be exactly satisfied, then a true complementary virtual work approxi- 
mate solution can be achieved. 

However, there may be problems for which it is too difficult to achieve satisfaction of 
all the constraints. Or it might be desired to generalize the formulation in order 
to achieve results which could be adapted to different sets of conditions. For what- 
ever the reason, in Equation (1), all or part of the boundary can be converted to 
boundary. The effect of this is to increase the number of displacement degrees of 
freedom on S^; that is, the magnitude of F in Equation (10) is increased. This modified 
problem can then be solved by complementary virtual work. The final stresses for the 
modified problem will be functions of the additional degrees of freedom, and these 
additional freedoms can finally be selected to impose some measure of satisfaction of 
the original requirements on the original surface. 
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SECTION 2 

DEVELOPMENT OF THE ASSUMED STRESS STATES 
2J INTRODUCTION 

In Section 1 it is shown that a complementary virtual work solution requires a 
stress state which satisfies equations of equilibrium in the domain and traction boundary 
conditions on the S^ portion of the boundary. For the finite element shown in Figure 1 , 
the exterior boundary is an S^ boundary because all displacements are specified to 
guarantee adjacent element continuity; therefore, there are no requirements on the 
stress state evaluated at the exterior boundary. 

The nature of the inner, or hole, boundary will depend upon the particular problem of 
interest. The inner boundary will be completely S^ for an open hole. If the hole is 
to be filled with other finite elements and complete adjacent element continuity is 
required, then the inner boundary will be completely in the complementary work solu- 
tion for the hole element. If the hole is to be filled with a fastener and if there is load 
transfer with separation between fastener and hole boundary, then part of the boundary 
might be completely and part will be completely S^. Or if the hole is to be filled 
with a fastener and if friction between fastener and hole boundary is ignored, then all 
or part of the boundary will be with respect to the radial direction (continuity of 
radial displacement between hole boundary and fastener) but with respect to tangential 
direction (zero shear stress on the hole boundary). 

The most general hole boundary conditions would be completely S^. Then any 
conditions could be satisfied by using the procedure discussed in Section 1 .7. 

However, all intended applications of this finite element will have zero shear stress on 
the hole boundary (either open hole or friction free fastener); and for these problems, a 
more accurate approximate solution should be obtained by requiring that the assumed stress 
states give exactly zero shear on the inner boundary. Therefore, this boundary condition 
is imposed on the solution which follows. 

Because the hole boundary is circular, the boundary condition can be most easily 
expressed through the use of polar coordinates. Therefore, in what follows, the domain 
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will be described by polar coordinates with origin at the center of the circular hole 
boundary (see Figure 1); and the stress state will be expressed in polar components. 

2.2 THE MOST GENERAL PERIODIC STRESS STATE 

In polar coordinates, the equations of equilibrium for plane stress and no body forces 
are Cl,p.66] 



ar 





2 

r 





(26) 


These equations are satisfied if the stress components are defined as follows in terms 
of a stress function 0: 


r 

r 


= i + 

r ar 


J_ 

2 

r 



r 


d 



(27) 


^r0 ar V 


Begin the development of the assumed stress state by recognizing that the stresses 

must be periodic in 0 in order to provide traction continuity across interior surfaces 

in the domain. Since any periodic function can be expanded in a Fourier series, it 

follows that the shear stress, T „ , can be written in the form 

r0 

00 

T^g(r,e) = f(r)+ 2 
n = 
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Ftxjtn this it follows that the stress function, 0, has the form 

00 

0(r , 0) =|^rF (r)J 9 + rG(0) + H(r) + ^ cos n 0 + B^(r) sin n 0j 

n “ 1 


Then T=9 (r F) + (periodic terms) 
® dr^ 


and since must be periodic, it follows that 
u 


(rF) = 0 


dr 


which gives F(r) - ^2 7 


Now r = - (—§- + G + G,0| + (periodic terms), from which it can be concluded that 

\d0^ ‘ / 


+ G(0) + C ^ 0 = Cg + ^ (a^ cos n 0 + b^ sin n 0) 

n = 1 " 


From this it can be shown that 


G(0) ~ 0 +C 2 + C^ cos 0 +C 2 sin 0 + ^a^ 0sin0- ■^b^ 0cos 0 


^ — ;r^ — (a cos n 9 + b sin n 0 ) 

^ 2 , n n 

n — 2 


Finally, it can be shown that the most general stress function which provides periodic 
stresses can be written as 
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(28) 


0(r,0) = H(r) + 0 ^ ^2 D2r0cos 0 

CO 

+ ^ cosn0 + B^(r) sin n sj 

n = 1 

Note that this stress function is valid for any material because the expression has 
been derived only from equilibrium considerations. 

2.3 RESULTANT FORCES ON THE CIRCULAR BOUNDARY 

Based on Equations (27) and (28), the stress components r and T „ are 

r ro 


_ 1 dH , 2 cos 0 2 sin 0 

T = - -j- + D„ D.J 

r r dr 2 r 3 r 


I 

n = 1 


n 

r dr 


2 

n 

r 


\r dr 


- jsin n 0 


(29) 


r0 


03 

= D, ^ n I (—A ) sin n 0 - -^ (— B ) cos n 0| 

12 ^ dr r n dr ' r n 

r _ , ^ 

n - 1 


(30) 


These stresses, acting on the hole boundary with radius of R, have a resultant force 
with components as follows (see Figure 2) 

.277 


- / (Tq sin 0 - cos 0) Rd 0 

2it 

F - ~ / (r sin 0 + r _ cos 0) Rd 0 
y J f r0 

•'0 


(31) 
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Substituting Equations (29) and (30) into Equation (31) gives the results 

F =-27TD_, F =2ffD, (32) 

X 2 y 3 

Note that the results shown in Equation (32) are independent of the magnitude of R. 

2.4 SATISFACTION OF TRACTION BOUNDARY CONDITIONS 

The traction boundary condition requires that t^q = 0 for every value of 0 at r = R. 

From Equation (30) it follows that 


at r = R, j (1 a J = 0 B ) = 0, n = l,2,... (33) 

dr r n dr r n 

The constant D, must be discarded; but the functions A (r) and B (r) are as yet 
I n n 

unspecified, subject only to the requirements of Equation (33). 

It can be shown that if Equation (33) is satisfied, then the resultant force components 
are still given by Equation (32). 

2.5 ORTHOTROPIC STRESS-STRAIN RELATIONS IN POLAR COORDINATES 

The finite element is constructed of a linearly elastic, orthotropic material with 
principal axes of orthotropy in the x, y coordinate directions. If initial stresses and 
strains and thermal strains are neglected, the stress-strain law in its most simple form is 
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with 


V V 

xy _ yx 

E E 

y X 


Note the order of subscripts in the definitions of the Poisson's ratios. 

Because the stresses and strains will be specified in polar coordinates/ it is useful to 
transform the stress-strain relations to polar coordinates. Begin with the transforma- 
tions of stress and strain, as follows [l , p.67]: 




^r0 


'Kyj 




xyl 


r0> 


(35) 


where 




2^ 
cos 0 

sin^0 


sin^0 

2 

cos 0 


L-2sin0cos0 2sin0cos0 


sin0cos0 

-sin 0COS0 

2 2 

cos 0-sin 0' 


(36) 


Then 




r0 


'r0 


r0; 


The desired material property matrix is / which has the following form; 

[Se ] = ^1 [''re] " S bU * S bU ^ =4 [<el 


(37) 


(38) 
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^1 8 Ie E E G 

\ X y y xy> 


C =1 (J-+-l+2^+ -L 

^2 2 I E E ^ E G 

\ X ■ y y xy, 


^3 2 Ie e 

\ X y 


c =I(± + -L+2-^--L 

4 8 Ie e e g 

\ X y y xy 


1 1 
1 1 

.0 0 


cos 20 


Se = ° 


0 -sin 20 

-cos 20 -sin 20 


l-sin20 -sin 20 


cos 40 

-cos 40 

-2sin40' 

-cos 40 

cos 40 

2sin40 

-2sin40 

2sin 40 

-4cos40 

.1 array is written in the form 


11 

^12 

^13 

21 

C22 

^23 

31 

^32 

^33 
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then 



= 

C ^ ^3 ^4 

C 22 

= 

C.J - Cg COS 20 + C 4 cos 40 

^33 

= 

r - 4 cos 40 

v -2 *-4 

^12 "^1 “ iS’S 

Si 

= 

^31 

= 

C ^2 ~“^3 20 - C 4 2 sin 40 

^32 

= 

C 22 ~ "Cg sin 20 + C 4 2 sin 40 


For possible use. Equation (39) can be inverted to give 


X 

1 


C,+C3+C4 


C] “ Cg +C^ 


Ik 

E 

y 

1 


xy 


■^1 2^2 "^^4 


C2-4C4 


Finally, note that for isotropic material, - 0. 


( 44 ) 


(45) 


2.6 COMPATIBILITY CONSIDERATIONS 

With an exact solution to any plane stress problem In polar coordinates, there is a 
necessary condition that the strains satisfy the compatibility equation 





+ r 


aVg) 

br^ 


S^(ry n) 

rd' 

br b0 


= 0 . 


(46) 
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With a complementary virtual work formulation, the strains will be derived from the 

assumed stresses by means of Equation (37); and there is no admissibility requirement 

that the strains should satisfy Equation (46). Indeed, the approximate satisfaction 

of Equation (46) is one of the consequences of the vanishing of the complementary 

virtual work. However, a better approximate solution should be obtained if the 

assumed stress state is compatible; so there is motivation to attempt to select the 

functions H(r), A (r), B (r) in Equation (28) so that the resulting stresses will be 
n n 

compatible. This was attempted but with no success. Therefore, in what follows, 
the stress components are not compatible for orthotropic material; and it is left up to 
the principle of complementary virtual work to achieve the "most compatible" stress 
state within the assumptions. 

The domain of the finite element is doubly connected, and it is known that Equation (46) 
is not sufficient to guarantee single valued displacements in a multiply-connected domain. 
Based on experience with isotropic material Cl, Sec. 43], it is anticipated that the 
displacements associated with the constants D 2 and Dg in Equation (28) will be multi- 
valued. These constants should not be discarded because they account for the non-zero 
resultant forces on the hole boundary, as shown in Equation (32). Therefore, additional 
terms should be added to the D2 and coefficients in order to enforce single-valued 
displacements. This will be demonstrated for the D 2 coefficient, as follows. Assume 
the stress function 


0(r/0) = D 2 r0 sin 9 + Arinr cos0 

Note that the additional term has the proper form to provide periodic stresses, as shown 
in Equation (28). Then from Equation (27) 


_ _ 2 cos 0 . A cos 0 

r - D., + A 

r 2 r r 

_ _ A COS0 

T. 

r0 r 
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with the strains derived from Equations (37), (43), and (44). Polar coordinate strain- 
displacement relations are [1 1 P»76] 


_ 5u 

S 3r 

_ u , Jl_ 

^0 r r S0 

_ 1 Su , V 

^r0 r "50 5r r 


( 47 ) 


Following the procedure outlined in [1, Sec. 31], the requirement for existence of 
displacement functions u(r, 0 ) and v(r, 0 ) can be written as 

Cf(0; C 3 , C^; A)] Inr - [A(4C^ - 2Cj + (4C^ - C^lsm 0 

+ g(C^, C^/ C 3 , C^; 03 / A) sin 30 + h(C^, 03 / C 3 , C^; 03 / A)sin50 

+ 11^ + F(0) + r^^ - G(r) = 0 ( 48 ) 

d0 

where F(6) and G(r) are functions which appear in the equations for displacement 
components u and v. 

Clearly Equation (48) can never be satisfied unless the function f(0; Cg, C^; 03 / A) = 0; 
and it can be shown that this is possible only for an isotropic material with €3 = 0 ^ = 0 . 
This is simply a manifestation of the fact that the stress state derivable from Equation (28) is 
not compatible for arbitrary values of the coefficients. 

If the Inr term In Equation (48) is disregarded, then the equation can be satisfied by 
requiring 


r^-G(r) = -K 


(49) 
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+ F(0) = K + CA(4C^ - 2 C 3 ) + D^(4C^ - C 2 )]sin 0- g sin 30- h sin 50 (50) 


Equation (50) will possess a multivalued solution proportional to 0cos0/ which would 
mean multivalued displacement components, unless the coefficient of sin 0 vanishes. 

Therefore, it is necessary to require 

1 

4Ci-C 

A = -D — 

^ ^2 4 C^ - 2C3 

_ " ^2 
Define ^ ^ 4 c - 2C 
I 3 


Conclude that the stress function 0(r,0) must include the following terms in order to 
eliminate multivalued displacements. 

0(r,0) = D 2 (r0sin 0 - ^rlnr cos 0) + D 2 (r 0 cos 0 +C rinr sin 0) (52) 

For isotropic material, C= ~~2~ / which agrees with Cl, Sec. 43]. 

2.7 FURTHER SPECIFICATION OF THE STRESS FUNCTION 

Equation (28), as modified by Equation (52) and = 0, includes arbitrary functions 

H(r), A (r), B (r) which are as-yet unspecified. Even though the resulting stresses 
n n 

cannot be compatible for orthotropic material, it appears reasonable to select these 
functions so that the resulting stresses would be compatible for isotropic material. 

Then the solution should be very good for the limiting case of isotropic material because 
the stresses will satisfy equilibrium and compatibility in the domain and traction 
boundary conditions on . Furthermore, if the orthotropic material is not too far 
removed from isotropic, in some vague and unspecified manner, a relatively few stress 
terms may provide a suitable answer. 

The isotropic stress function is shown in [1 , p.l33] which gives the following expressions 
for the functions: 
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(53) 


H(r) = Qq Inr + bg 

A, (r)=b,r^+aJ 1 

B, (r) = d,r^ + c, i 


(54) 


A (r) - b r + a r 

n n n 

B (r) = d r + c r 

n n n 


, , -n+2 -n 

+ b r + a r , 
n n 


I 


+ d r 
n 


-n+2 


+ c 

n 


-n 

r 


. n = 2, 3, . . . 
n = 2, 3, . . . 


(55) 


However/ there is a flaw - possibly a fatal flaw - in the use of Equation (55) for 

A and B . It will be seen later that in order to obtain a sufficient number of stress 
n n 

coefficients it might be necessary to include terms up through A^^/ B^^ or possibly 

1 6 

even more. This would mean that the radial dependence would include terms from r 
-12 

down to r . This extreme range of exponents is likely to lead to a very ill- 

conditioned array CH], which is derived as shown in Equation (14) by integrating over 

the domain of the finite element. For this reason, Equations(55) are discarded, and 

the functions A (r), B (r) are assumed as follows: 
n n 

2 I II 

A2(r) = 02r 

(56) 

» 2 ( r ) ' '*2 =2 T 

r 


r 

= =3''^ ■^3 7 =3 -3 

r 


(57) 


26 



N 


M 


A(r)= y C?r'+ y C" -L, 


= 0 
N 


j - 1 - r- 

M 


B„W =2 S|' r' + y s" -L , 

n • _n • - I J i 

1=0 J = 1 r J 


n = 4, 5, ... 


n= A, 5, . . , 


( 58 ) 


In Equation (56), the r term has been neglected since^this term vanishes for the case 

of an infinite isotropic plate under uniform stress at infinity Cl , p. 91 ]. In Equation (57), 

5 

the r term is neglected because it does not appear in the case of an infinite isotropic 

plate with linear variation of stress at infinity C2, p. 470]. Finally, Equation (55) has 

been discarded completely for n ^ 4 and replaced by Equation (58), which expresses 

A (r), B (r) as polynomials in r. The functions are written in a form to allow for differ- 
n n 

ing numbers of terms with positive and negative exponents. This flexibility in specifying 

functions A (r) and B (r) will permit the selection of a sufficient number of stress coef- 
n n 

ficients without the extreme range of exponents which naturally appears with fully 
compatible stresses. Note that the resulting stress state will now be incompatible 
even for isotropic material. 


This numerical difficulty with compatible stresses is obvious in Equation (55) which,of 
course, is valid only for isotropic material. It is suspected, however, that the same 
problem will arise if orthotropic compatible stresses are used for the assumed stress state. 
This suggests that a fully compatible solution for orthotropic material might never be 
achieved due to appearance of ill-conditioned matrices in the formulation. 


2.8 CONSIDERATIONS OF SYMMETRY 

If two coordinate axes, say x and y, are established in a plane, there are four possible 
symmetry conditions which can be established relative to those axes. These conditions 
can be defined by describing the following four functions which possess the symmetry 
properties: 


SS(x,y) = SS(x,-y) = SS(-x, y) = SS(-x, -y) 
SA(x,y) = SA(x,-y) =-SA(-x,y) = -SA(-x,-y) 
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AS(x,y) = -AS(x,-y) = AS(-x,y) - -AS(-x,-y) 


AA(x,y) = -AA(x,-y) = -AA(-x,y) = AA(-x,-y) 

The function SS(x,y) is symmetric about both x and y axes; it will be called a 
symmetric-symmetric (SS) function. Function SA(x,y) is symmetric about the x-axis 
and antisymmetric about the y-axis; it is a symmetric-antisymmetric (SA) function. 
Likewise, function AS(x,y) is an AS function; and AA(x,y) is an AA function. 

With the four symmetry conditions established, it can be shown that any function f(x,y) 
can be written as the sum of four functions, each of which possesses one of the four 
symmetry conditions. This is demonstrated by first writing the identity 

f(x,y) = ^[f(x,y) +f(x,-y) +f(-x,y) +f(-x,-y)] 

+ j Cf(x,y) + f(x,-y) -f(-x,y) -f(-x,-y)] 

+ ^ C f(x,y) - f(x,-y) + f(-x,y) - f(-x,-y)] 

+ j Cf(x,y) - f(x,-y) - f("X,y) +f(-x,-y)] 


Define SS(x,y) = j Cf(x,y) +f(x,-y) +f(-x,y) +f(-x,-y)] 
SA(x,y) = ^ [f(x,y) +f(x,-y) - f(-x,y) - f(-x,-y)] 
AS(x,y) = Cf(x,y) - f(x,-y) +f(-x,y) - f(-x,-y)] 
AA(x,y) = j ljF(x,y) - f(x,-y) - f(-x,y) + f(-x,-y)] 
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It can be confirmed that these functions do indeed possess the indicated symmetry 
properties and, as stated, 

f(x,y) = SS(x,y) + SA(x,y) + AS(x,y) + AA(x,y) 


What has been described above for functions of x and y is also applicable to stress 
and displacement states. That is, any arbitrary stress and displacement state can be 
represented as the sum of SS, SA, AS, and AA states. These four states are defined 
in Figure 3. Note that a stress state is symmetric about a line if the stress states at 
two symmetrically located points are mirror images of each other; the stress state is 
antisymmetric about a line if the stress states at two symmetrically located points are 
negative mirror images of each other. Similar definitions apply to displacements. 

These physical definitions of symmetry and antisymmetry for stresses and displacements 
mean that individual stress and displacement components might actually have different 
characters when signs are established according to usual conventions. For example, 
in the SS state of Figure 3, the stress T is an SS function while stress T « is an AA 
function, displacement is an AS function, and displacement is an SA function. 

2.9 FINAL FORM FOR THE STRESS FUNCTION AND STRESS COMPONENTS 

The radial coordinate, r, is first normalized or nondimensionalized in terms of the 
radius of the circular hole, R, giving 


r 

P= D 


R * 


Then, if the stress function is written as 0(p , 0), the stresses will be given by 


r2^ = 1 + _L 

^ P 302 


R^r = ^ 
R Tg ^ 

3p 


( 59 ) 




r0 3p Vp 90 
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Figure 3a . Symmetric-Symmelric State 
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Figure 3b. Symmetric - Antisymmetric State 
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Figure 3c, Antisymmetric-Symmetric State 
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Figure 3d. Antisymmetric-Antisymmetric State 
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and the equilibrium equations are 






( 60 ) 


The stress component T « is first calculated from the stress function given by 

r9 

Equations (28), (52), (53), (54), (55) with = 0. The arbitrary constants are 
then related as necessary to satisfy the traction boundary requirement that T^q~ 0 
at p = 1 . The final stress function is then divided into its four components according 
to the four possible symmetry conditions of the stress state, with the final result shown in 
Equations (62) and (63). The three stress components are derived from Equation (59) 
and recorded in Equations (64) , (66) , and (68) . It can be confirmed that these stresses 
do indeed satisfy the equilibrium equations. Equations (60), and provide at p= 1 . 

In Equations (62), (64), (66), and (68), the SS, SA, AS, AA symmetry conditions 

th 

are numbered 1, 2, 3, 4 respectively. The coefficient B. . denotes the i coefficient 
th 

for the j symmetry condition. The quantities Nl,N2,N3,N4 denote the upper 
limits on the harmonic expansion for each symmetry condition. Then for each harmonic 
above the third, the quantities II, 12, 13, 14 and Jl, J2, J3, J4 denote the upper 
limits on the summations of terms with radial coordinate dependence. These quantities 
could be functions of harmonic number; but in what follows, they are assumed constant. 


For harmonics beyond the third, there is obviously a great deal of arbitrariness in the 
number of terms included in the assumed stress state. However, it is possible to number 
the coefficients in such a way as to allow for the arbitrariness. 


Define Kin = 4 + (^)(H + Jl) , K3n = 4 + (^^^)(I3 + J3) 


K2n = 4 + (^)(I2 + J2) , K4n = 2 + (^)(I4 + J4) 


( 61 ) 
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With the definitions of Equation (61), it con be confirmed that the coefficient subscripts 
shown in Equation (62) will be properly sequenced for arbitrary values of N,I,J. 


where 


0(p,e) = B^^(M + B2^(^P^ 
^ ^31 I 


(1 + p^) + B^^ ^ ^ ^ 


SS 


N1 

^ 2 


II ji 

2 ,F(p;i)+ 2 B,, 


n^,6li=l --Kin,! — .^^>Kln.Il,l 


G(p;j) 


cos n0 


+ Bi 2 |p0 sin 0 - C (p Inp- ^ p^ cos 0^ ^ ^ ® 

^[^32 I ^ ^ + B^2 ^ 2p^j cos 30 


12 


N2 

2 

n =5,7 Li=l 


J2 


2 ®i+K2n,2 ■*■ 2 VK2n+I2,2^^'j^ 


j=l 


cosn 0 


SA 


B^ 2 jp 0 cos 0 + C (p Inp - P^ sin 0^+ B22 (p ^ + P^ sin 0 
[®33 ^ ^ '^^43 'll ^ 


N3 


^2 


13 


J3 


n=5 7|_^i ®i+K3n,3 VK3n+I3,3 


sin n 0 


■ [^Bi 4 ^ (1 + P^) + ^ ^24 (p 3p^)j sin 2 0 


N4 

^ 2 

n - 4,6 


14 


2 


J4 


“i+K4n,4 VK4n+I4,4 


F(p;i) = p' + i-l 
G(p;j) = p ^ - j + 1 


AS 


AA 

sinn0 (62) 


(63) 
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R^r^(p, 0 ) = B^^(p"^) +B 2 ^ 0 ) 


“|b3i( 1 +2p + B^^(l +p ^)jcos 20 


N1 

^ 2 

n=4,6 


II 


J1 


i?i ^'+Kln,l 2 ^ VkIh+II,! 


cos r 


+j ^12 ^ ^ ■*■ B22(P"P ^ j cos 0 

- |b 22 (3p+ 5p ^ + B^2 P ^)j ‘=°s 30 


N2 

^ 2 

n=5,7 


12 


J2 


Vk 2 o ,2 ^j+K 2 n+I 2,2 


cos 


+ |b^31^-2p ^ + C (p ^ - p)j + B 23 (p-p ^| sin 0 
- l^Bgg (3p + 5p ^ + B^gCp + ft ^)j SIP 30 


N3 


n=5,7 Li = l 


13 J3 

S+K3n,3'r'^'''*" ' “j+K3n+I3,3'"r 


^ W 

2 3 ; 4 -k"?„ T F^(pM/n) + 2 JP'j'^) 


sin n 


-[' 


8,4(1 +2p'^)+B24 


(1 + p '*) si 


sin 20 


N4 

^ 2 

n-4,6 |_ 


14 J4 

®i+K4n,4 ® +X4n+I4,4 


sin n 


where 


F/p; i, n) = (i - n^)p'"^ - n^(i - l)p'^ 
G|.(p;j/n) = -(j+n^) p +n^(j+l) p ^ 


SS 

0 


SA 
n 0 


AS 

0 

AA 

0 (64) 

(65) 
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R%(P/0) = B^^(-p'^)+B2i (1) 


+ [®3] ®41 P “s 20 


SS 


N1 

^ 2 


II 


ji 


2 1 Fo(p?0 + 2 


6'x .-1 I+KlnJ .e, >Kln+n,l ^0 

n— 4/0 Li “ ' J — ’ 

+{b ^2 [*“ ^2p- p + B22(3p + P cos 0 
+ [B32 (3P + P ^ + B^2 (P + P cos 30 


Go(p;j) 


cos n 0 


SA 


N2 


n=5/7 


12 . J2 

;?i ^'+K2n,2 *"0^^''^ ■*■ VK2n+I2/2 


cos n 0 


+{b^ 3[- C (3p-p“^] + B23OP+ p'^} sin 0 
+ [B33 (3p+ p ^ + B^^Cp + p sin 30 


AS 


N3 r 13 J3 

2 2 B;+|<3„^3 Fg(p;i) + 2 

n-O// L * 


B 


j+K3n+I3,3 '^0 


GJp;j) 


sin n 0 


+ [b ^4 (D+B 24 (l+p“'^)]sin 20 


N4 

^ 2 

n-4/6 


14 

2 


J4 


Fa(p;i) + 2 B^^ > i Gp(p;j) 


Li = l i+K4n,4‘0'^'’^ t^i j+K4n+I4/4^0 


AA 

sin n 0 (66) 


where Fg(p; i) - i(i - 1) p 


i-2 




(67) 
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where 


“■ [®31 ■ P ®41 ■ P "*)] sin 20 


N1 


n=4,6 


II 

,t, ®i+Kln,l .2 Vln+11,1 


J1 


ss 


sin n 6 


+ {b^ 2 [^ (p-p“b] + B 22 (P-p'^}sin 0 
+ [632 (3p- 3p"^ + B^ 2 (P-p'^)] sin 30 


SA 


N2 

n=5,7 


12 


J2 


,1, ®i+K2n,2 V2n+E,2 


sin n 0 


+{b^ 3 [c (p~p B 23 (P-P ^j-cos 0 
- l^Bgg (3p- 3p ^ + B^g(p “ P ^)] COS 30 


AS 


N3 

- 2 
n=5,7 


13 J3 

.|/i+K3n,3^e<^’''">^ 2, ®J+K3n+13,3 ® 


■ [®14 ®24^’ ■ 

14 J4 

®i+K4n,4 2 Bj+|< 4 n+i 4^4 Gj.Q(p;j/n)J cos n 0 


N4 

- 2 
n=4,6 


AA 


(68) 


F^g(p;i /n) - n(i - l)(p‘ ^ - p ^ 
G^g(p;j,n) = n(j + 1)[p"^ - 


(69) 
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SECTION 3 


SPECIFIED BOUNDARY DISPLACEMENTS 
3.1 EXTERIOR BOUNDARY 

As mentioned earlier, the exterior boundary is completely S^ type, with displacements 
specified so as to maintain continuity with the adjacent elements. In the work which 
follows, it is assumed that each adjacent element contains only two nodes on the boundary 
common with the hole element, with linear displacement between nodes; see Figure 4a. 
Therefore, the displacement on the exterior boundary of the hole element will be 
specified as linear between nodes. If a nondimensional coordinate ^ is established 
such that ^ = 0 at node i and 4 = 1 at node j, then both displacement components 
between i and j will have the form 

v(|)=v.(l-5)+v,(«) (70) 

where v. denotes the value of displacement v at node i. These specified displacements 
obviously satisfy the rigid body mode requirement of Equation (21) because the most 
general rigid body motion of an initially straight segment between nodes will result in 
linear displacements. 

If the adjacent elements have a more complicated boundary displacement, this 
complementary virtual work formulation can very easily be modified. For example, 
suppose the adjacent elements contained three nodes per boundary common with the 
hole element and suppose the displacement is quadratic along that boundary; see 
Figure 4b. Then the displacement components on the hole external boundary between 
nodes i,j,k will have the form 

v(5) = v.5(i^-5)+Vj(l-|^)+V|^i(5^ + 5) . (71) 

If this portion of the boundary receives a rigid body displacement, then it will 
follow that 

v^4(v,+v,) (72) 
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Then if Equation (72) is substituted into Equation (71), the displacement at every 
point between nodes i and k is given by 

v(^) = Vi^(l +?) 

This linear displacement function shows that Equation (71) contains rigid body motions 
and does satisfy Equation (21). 


3.2 INTERIOR BOUNDARY 

Although the circular hole boundary is with respect to shear stress, it is with 
respect to radial displacement. Therefore, it is necessary to specify the radial dis- 
placement, u^, on the boundary. The generalized displacements will be the values 
of displacement at discrete points (the nodes) with displacements between nodes 
given by specified functions of position on the boundary. 


Consider a segment of the boundary between two nodes; see Figure 5. 


Let 0. = initial angle for the segment 

0^ = final angle for the segment 

Begin with the assumption 

u^(0) = u.g.(0) + Uj.g^(0) 
subject to the requirements that 

g;(0.) = 1 / 0i^®f^ ~ ° 

gf(0i) = o / gf(0f) = ^ 

which guarantee that Equation (73) does furnish the proper nodal displacements. 


(73) 


( 74 ) 
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a) Linear Edge Displacements 



Figure 4. Displacement Continuity Between Adjacent Elements 



Figure 5. Segment of Hole Boundary Between Two Nodes 
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Now consider a rigid body motion of the hole boundary/ with components u , v 
in the x, y directions. The radial component of this rigid body motion will be 

u^(6) = u cos 0 + V sin 0 (75) 

which means that the nodal values of rigid body displacement are 

u. = u cos 0. + V sin 0. 

I I I 

(76) 

u^ - u cos 0^ + V sin 0^ 

If these nodal values are substituted into Equation (73) the displacement between nodes 
is given by 


or 


u^(0) - (u cos 0. +v sin 0.) g.(0) + (u cos 0|. +v sin 0|.) gj;(0) 

Ur(®) ~ “fsi (®^ 9f(®) f^i (®^ ®i ^f(®^ ®f] 


If Equation (73) is to contain rigid body modes, it is necessary that Equation (77) have 
the same form as Equation (75) which represents rigid body motion between nodes. 

This leads finally to the requirement that 


cos 0. 

cos 0^ 


CD 

O) 


COS 0 1 

sin 0. 
1 

sin 0J. 


9f(e) 


sin 0 I 


which gives the following unique expressions for g.(0) and g^(0): 

sin (0^-0) sin (0-0.) 

9i^®^ " sin (0^-0.) ' " sin (0^-0.) 


( 78 ) 


It can be seen that Equations(78) satisfy Equations (74). Therefore, the final form of 
the specified radial displacement between nodes is as follows: 
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(79) 


sin 

u (0 ) = U. : /q 5-t 

r I sin (0^- 0.) 


+ u 


sin (0- 0.) 
f sin(0^- 0.) 


This assumption provides proper nodal values and does contain rigid bod/ motions, 
which guarantees that Equation (21) is satisfied on the hole boundary. 


If the boundary nodes are given the same radial displacement, it would be desirable 
for the displacement u^ to be uniform between nodes. In order to see if Equation (79) 
contains uniform displacement, set u.= Uj. = u and get 


" " sin(e’-ep +»!n(e-e.)] 


(80) 


Unfortunately, Equation (80) does not reduce to uniform displacement, which means 
that Equation (79) does not contain uniform displacement. However, if (0|.-0.)«1, 
it can be shown that Equation (80) gives the approximate result 

u (9) = u 
r 

Therefore, Equation (79) will give essentially uniform displacement as the arc 
length between nodes becomes sufficiently small. 


The following displacement function does include uniform displacement for any arc 
length: 


° 2sin('^-e.) ‘ V] 

2sin(9y-ft) t'" <®f • • ®i >] 

However, Equation (81) does not include rigid body motions. 

in conclusion, it is impossible to devise a displacement function which contains both 
rigid body motion and uniform displacement. It is much more useful to include rigid 
body motion, so Equation (79) will be used in the following work. 
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SECTION 4 

DEVELOPMENT OF THE ELEMENT STIFFNESS MATRIX 


4.1 STRESS STATE DUE TO IMPOSED BOUNDARY DISPLACEMENTS 

Begin with the expression for complementary virtual work, as given in Equation (7) 
and repeated below: 


CVW = 


= [Srj {c } dA - y |_6tJ [u 3 ds 


In polar coordinates (Figure 1), the differential area element can be expressed as 
follows in terms of the nondimensional P = "^ 


dA = rdrd0 = R p dp d9 


and the differential boundary length can be written as 


ds = Rf(6)de 

where the function f(0) depends upon location on the boundary. 

0 2 
The real strains, {R^c}, ore expressed in terms of real stresses, [R t3, through the 

stress-strain relations in polar coordinates, as given in Equations (43), (44), and (39): 


tR^t} = CC^3] 


Then the CVW can be written as 


CVW 


= 1-J Lr^6tJ Cc^] {r^t 3 pdpde - ^ i?3f(0)de 


NSC 


Assume [R 6 t 3 - y CAJ(p,e)]Co^] 

3=1 



2 ^ i 

CR 6T} = ^ ItHQ)! [a^} 
j=l 


, NSC 

[rM= 2 

j = l 


CA^(P, 9)] [fij} 


NSC 

|u}= 2 CC^(0)][q^3 

2 = ^ 


These assumptions are consistent with what was shown earlier in Equations (8), (10), 
(11), and (12), with expanded notation to recognize that each stress, surface traction, 
and displacement vector can be written as the sum of four symmetry condition vectors. 
In Equations (85) through (88), the superscript j denotes symmetry condition number; 
and NSC denotes the number of symmetry conditions to be considered in a particular 
problem'. In general, NSC will be equal to four; however, if the problem is known 
to contain only certain symmetry conditions, then NSC can take on the proper and 
limited set of values. 


If Equations (85) through (88) are substituted into Equation (84), the result is 
NSC NSC , . r . T 


NSC NSC I ,|/r. r . \ n r 


} 


i = 1 k = 1 



R S 


[r^] cr]f(0)d0 


it) 


Now, the domain and boundary are symmetric about both the x and y axes. There- 
fore, a boundary integration of the product of symmetric (antisymmetric) tractions and 
antisymmetric (symmetric) displacements will vanish. Likewise, a domain integration 
of the product of symmetric (antisymmetric) stresses and antisymmetric (symmetric) 
strains will vanish. Also, note that because the x, y axes are material'axes for 
homogeneous and orthotropic material, it follows that symmetric (antisymmetric) 


stresses will result in symmetric (antisymmetric) strains. Therefore, because of the 
symmetry properties of domain, boundary, and integrands, it follows that 


JCA''3pdpde = 0, if 


Cr]f(e)d0 = o, 


if 


f lA^l iC dp 66 = 4 [ CC JCA^lpdpde 

JA ./Quadrant 

L Wim 66 = 4 1 Crj3‘‘’ccj]f(e)d0 

./Quadrant 


Define 


= / [AJ] [Cg 
J Quadrant 


] Ca j] p dp d0 


CGj] = / crj]'*'cc^]f(0)d0 

•/Quadrant 


After substituting Equations (90) through (93) into Equation (89), the CVW can be 


written as 


CVW = 4l ^ [cjjjl [hJ: {b 1} - [Gi] tq^l 


If CVW = 0 for every [a^J , conclude that 


{ - CQj] {qj} = CO}, j = 1, 2, 3, 4 
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M=RCHj]“’ CG^] [qj] , j=l,2, 3, 4 

Define [60^] = CH^]”’ [G^] , j = 1 , 2, 3, 4 

Therefore, [B^] = R [BGi^] {q^} , j = 1,2, 3, 4 

The final expression for the stress state is found by substituting Equation (95) into 
Equation (87) and solving for [t] to give 


NSC 


{r(p,9)} = i f CaJ(p,9)]CBqJK5J} 

j = I 


(94) 

(95) 


(99) 


Equation (96) will give the stress state at any point in the domain due to imposed 
boundary nodal displacements. Note that the nodal displacements are expressed in 
terms of the four symmetry condition components. 


4.2 DEVELOPMENT OF STIFFNESS- MATRIX 

In the anticipated applications of the hole finite element, the total or overall 
structural problem will be solved through a displacement formulation which can be 
most generally based on the principle of virtual work. Therefore, the hole element 
solution of Equations (94), (95), and (96), based on complementary virtual work, 
must now be placed into a form which is compatible with the formulation for the 
remainder of the structure. This is done by developing a stiffness matrix for the 
hole element which can then be handled in a routine fashion in the displacement 
based finite element programs. 


Begin with the expression for the internal virtual work done by an equilibrium real 
stress field moving through compatible virtual strains: 



( 97 ) 
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For use in a displacement formulation, the real stresses {f} must be expressed in 
terms of real nodal displacements; and the virtual strains {6c} must be expressed 
in terms of virtual nodal displacements. But this relationship between nodal dis- 
placements and stresses is precisely what is given, to the best approximation possible, 
by Equations (94), (95), and (96). Therefore, Equation (97) can be rewritten in 
terms of displacements by first introducing virtual strain-virtual stress relations in 
the form of Equation (34) and then substituting the stress-displacement relations of 
Equation (96): 

6W.=-~ / LR^6tJ [C .] CR^r}pdpd6 


aw. = - 
' 



CA^] Cc^g] CA'"]dA 


[BQ*"] Cq'^} 


(98) 


After accounting for the symmetry properties of domain and integrand, and after 
introducing Equation (92), it follows that 


NSC 


6W, = - 2 M 
i=1 


4tCBQj] CH^lCBOj] 


{q^3 


Define CKQ^]=t[BQj] CHj]CBQj] 

NSC 

Finally, 6W. = - ^ [aq^J UCKQ^] 

i= 1 




(99) 


( 100 ) 

( 101 ) 


After comparing Equation (101) with the result which would have been obtained 

through a virtual work development of 6W., it is clear that the array 4[KQ^] is 

^ th 

the stiffness matrix associated with displacements in the j symmetry condition. 
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More precisely. Equation (100) defines the complementary virtual work approximation 

th 

to one-fourth of the stiffness matrix associated with boundary displacements in the j 
symmetry condition. 

4.3 THREE TYPES OF HOLE ELEMENTS 

The hole finite element will generally be employed as shown in Figure 6; that is, the 
total square domain will be included In the finite element model of the structure. 
However, if the problem has known symmetries in structure and loading, then it is 
possible to achieve significant economies in both formulation and solution by working 
with only a portion of the structure. 

For systems which are known to possess only one of the four symmetry conditions, 
the simplifications begin with Equations (85) through (88) in which the summations 

over all symmetry conditions are eliminated; all quantities are expressed immediately 
in terms of only the one known symmetry state. Further simplification is achieved in 
the finite element model of the structure, as shown in Figure 7 for a typical case which 
is SS in both structure and loading. It is necessary to consider only one quadrant of 
the structure, which means that it is necessary to consider only one quadrant of the 
hole finite element when developing the stiffness matrix. 

Systems with symmetry or antisymmetry about only one axis can be considered to be 
the sum of the appropriate two symmetry conditions, which again simplifies Equations 
(85) through (88). Figure 8 shows a case with symmetry about the y-axis, which Is 
therefore a combination of SS and AS. The model can be simplified to only one half 
of the structure, which means that only one-half of the hole finite element Is necessary. 

Evidently there are three different types of hole elements which could be incorporated 
In a finite element model - the complete or total element in Figure 6, the half element 
in Figure 8, and the quarter element in Figure 7. With appropriate modifications to 
the symmetry condition summation. Equation (96) gives the stresses for each type of 
element; no adjustments are necessary. However, the stiffness matrix expressions in 
Equations (100) and (101) must be modified for the quarter and half elements. 


49 




50 




t t t 



I 1 I 

II Y 


Figure 8. Problem with Symmetry About One Axis, Permitting 
Use of the Half Hole Element 
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In Equation (98) it was assumed that the domain integration was taken over the 
complete hole element. Then 



)dA = 4 


i 


( ) dA 

Quadrant 


which accounts for the multiplicative factor of 4 in Equation (99). Obviously, if 
the structural problem is to be solved with the half element in Figure 8, then 



( )dA = 2 


i 


( )dA 
Quadrant 


r 


and with a quarter element, the physical domain is equal to the quadrant used for 
integration. 


It is now possible to give a physical interpretation to [KGr] in Equation (100). 

Clearly, [KQ^] is the quadrant element stiffness matrix associated with displace' 
th 

ments in the j symmetry condition. 


These results concerning stiffness matrices can be summarized in terms of [KQ^] 
of Equation (100), as follows: 


for a quadrant element with j 


th 


symmetry condition. 


CKQj] = CKQ^] (102) 

6W. = - [_6q^ CKQj] {q j} (103) 

for a half element with superposition of two appropriate symmetry conditions, 

CKH^] = 2[KQ'^], for two values of j (104) 

6W. = l_6qjJ CKH^] (q j] (105) 
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fora total element, requiring all symmetry conditions 


CKT^ = 4CKQJ] , j= 1, 2, 3, 4 

6W. =- ^ L^qjJ [KTj][qJ] 
j = l 

For possible future use. Equations (103), (105), and (107) can be written as 
6W. = [Sq^ [MCKQj]] { q j} 

where the summation is over appropriate values and M = 1 , 2, or 4 depending upon 
element type. 

4,4 TRANSFQR/V^ATION OF DISPLACEMENTS 

If the hole finite element is to be combined with other elements and a problem solved 
with a conventional displacement based program, then it is necessary that the nodal 
displacements for adjacent elements must be identical at nodes where continuity is 
required. This nodal continuity is most easily achieved if the nodal displacements 
for the adjacent elements are expressed in the same coodinate system and have the 
same physical meaning. In the following discussion, it is assumed that adjacent 
element displacements are in the same coordinate system (either x, y system or radial) 
so that it Is only necessary to consider the physical meaning of the measures of dis- 
placement. The two possible meanings are described as either the symmetry condition 
displacements or the general state displacements. 

Now, for problems with a single symmetry condition which are to be solved with the 
quadrant element, there is an identify transformation between symmetry condition 
displacements and the general state displacements. That is, if the system is, say, SS 
in character, then the general displacements will have zero contribution from SA, AS, 
and AA; and the general displacement will be identical to the symmetry condition 


(106) 

(107) 

(108) 
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displacement. In this case, the quadrant element stiffness matrix is given by 
Equation (100); no further modification is necessary, and the hole element quadrant 
stiffness can be immediately merged with adjacent element stiffness matrices. Then 
the problem can be solved with standard displacement based finite element programs. 

However, for the half element and total element, the displacements in Equations (105) 
and (107) are symmetry condition displacements; and these are not the same as general 
state displacements which represent the superposition of two or four symmetry condition 
displacements. Therefore, the symmetry condition displacements will be transformed 
to general displacements before attempting to merge stiffness matrices. The details of 
the transformations depend upon the number of boundary nodes, but the basic features 
will be presented for a simple case of each type of element. 

First consider the total element shown in Figure 9. Note that the nodal numbers begin 
at the first quadrant comer node and proceed counterclockwise around the element 
until all nodes on the outer boundary have been numbered. The next node number is 
assigned to the circular boundary node at 9 = 0®, and numbering continues counter- 
clockwise. Because the external boundary is entirely S^, there are two displacement 
degrees of freedom at each exterior boundary node. However, there is only the one 
radial degree of freedom at each interior node. On the exterior boundary, positive 
displacements are in the positive coordinate directions. On the interior boundary, 
positive displacements are outward. These sign conventions and conventions for nodal 
and displacement numbering are followed for all arrangements of nodes. 

Figures 10a, 11a, 12a, and 13a show the nodal displacement patterns associated with 
SS, SA, AS, and AA displacement fields, respectively. Note that the independent 
nodal displacements are selected to be in the first quadrant, with all other nodal dis- 
placements expressed in terms of the independent values. Also note that the quadrant 
nodal displacements are numbered beginning at 0 = 0® on the external boundary and 
returning to 0 = 0® on the inner hole boundary. This convention is followed for all 
arrangements of nodes. 
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a) SS Displacement Field for Total Element 



Figure 10, Symmetric-Symmetric Displacements for Total Element 





Figure 11. S/mmetric-Antisymmetric Displacements for Total Element 
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a) AS Displacement Field for Total Element 



Figure 12. Antisymmetric-Symmetric Displacements for Total Element 





Figures lOb, 11b, 12b, and 13b give the relationships between the independent 
nodal displacements and the general nodal displacements of Figure 9. As a check on 
these relationships, it can be confirmed that the summation of Figures 10b, 11b, 12b, 
and 13b, extended to the other three quadrants, will indeed give Figure 9. For 
example, consider the x-direction displacement of the first quadrant comer node. 

Then 

(qj) ss + (I 2 ) SA ^ AS ^ <‘l2> AA " 3 ("r V “9 * “ 1 3> ^ 

i(u, + Uj+ u,+ u,3> (u, - Uj -ujj) + ^(u, + Uj-u,- u,3> = u, 

as required. As one final example, consider the /-direction displacement of the exterior 
node at 0=180°. Then 

AS ■ <‘’l> AA “ 5 ^“8^ “16> • 5 “ “8 

as required. 

A half-element with symmetry about the y axis is shown in Figure 14. Node numbering 
begins at the lower left exterior comer, proceeds around the outer boundary, and returns 
to the inner boundary at <i>= 0°. There are two degrees of freedom at each exterior 
node except the nodes on the symmetry axis, which have only one displacement. 

Positive displacements are as shown in Figure 14 and described for the total element. 

The general state is a superposition of SS and AS due to assumption of symmetry about 
the y-axis; see Figures 15 and 16. Again, the independent displacements are selected 
to be in the first quadrant, with the same numbering scheme described earlier. As an 
example of the superposition, consider the radial displacement of the inner node at 
0 =45°. 

^VsS " ^‘^6^AS ^ 2 ^"l0‘^“l2^ “ '^““l0'*'“l2^ ‘^lO 

as required. 
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Figure 14. Node and Displacement Numbering for Half Element 
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b) SS DecomposiHon of General Displacements 


Figure 15. Symmetric-Symmetric Displacements for Half Element 




a) AS Displacements for Half Element b) AS Decomposition of General Displacements 

Figure 16. Antisymmetric-Symmetric Displacements for Half Element 


It is now clear that for each symmetry condition/ the general state nodal displacement/ 

Cu}/ and the symmetry condition nodal displacements/ {q^}/ are related by 

Cq^} = CAJ]{u) (109) 

where the transformation array depends upon the type of hole element (total or 
half) and the details of the nodal arrangement. Recall that for a quadrant element, 

CA'^] - CI]/ the identify matrix. 

4.5 TRANSFORMATION OF STIFFNESS MATRIX AND STRESS MATRIX 

The necessary transformation of the stiffness matrix is established by substituting 
Equation (109) into Equation (108) in order to express the internal work in general 
state displacements. 



6W. = ^ L® |mCkq^]J CaJ] [ u } 


Define 

T 

im = [AJ] 

[KQi]CAj] 

( 110 ) 

Then 

6w. = - [e uj 

tu} 


Define 

CK] = M ^ [K^] 

( 111 ) 


and it is clear that Equation (111) supplies the proper hole element stiffness matrix 
associated with general state displacements and this stiffness matrix is suitable for 
standard treatment by displacement based programs. 

Equation (96) gives the stress state in terms of symmetry condition displacements, 

{q^}. This equation can also be transformed to general state displacements as follows; 

Cr}=^X CAj]CBQj]CA^]{u} (112) 
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A more efficienf evaluation of the stiffness matrix is possible by substituting 
Equations (94) and (100) into Equation (110): 

.T .T . -T . .-1 

Ck^] = t Ca^] Cg^] Ch^] Ch^] Ch^] Cg^] [a^] 

Define CG^] =[G^]CA^] 

T . -1 

Then, CKj] = t[Gj] [H^] CG^] 

where it is recognized that array CH^], as defined by Equation (92), is symmetric. 
Equation (1 12) can also be reformulated as follows: 

-1 . . -1 . 

Define CS'^] = CH-^] CG^] CA-^1 = CH^] CG^] 

Then, {r(p, 0)} = ^ C (p, 0)] [S^] [ u } 

4.6 COMMENTS ON PROCEDURE 

The essential feature of the development presented in this section is the realization 
that symmetry conditions exist for stress, tractions, and displacements and then 
exploiting this feature to the fullest. This permitted the derivation of stiffness and 
stress matrices for the total hole element with the work integrations limited to just 
one quadrant of the element. This is a significant savings of computational effort 
compared to those formulations which integrate over the total element. 

When solving problems with the quadrant or half elements shown in Figures 7 and 8, 
the internal stresses are guaranteed to be continuous across the symmetry boundaries. 
This is true because the solution is in reality based on a total element which has 
correct symmetry conditions imposed; and the stresses are defined continuously in 


(113) 

(114) 


(115) 

(116) 


65 



the total domain. This can be contrasted with an alternative approach which begins 
with a true quarter element with no reference to the total domain. If the comple- 
mentary virtual work procedure is used to develop the quarter element stiffness matrix, 
and if four of these quadrants are then combined properly to model a total element, 
there will be traction discontinuities across those boundaries which are interior to the 
total element; and these discontinuities could cause troubles when attempting to 
evaluate stresses in the interior of the total element. 

However, it should be noted that the imposition of additional specified displacements 
along interior lines will tend to stiffen the total finite element when compared to an 
element with no specified interior displacements. The implications of this with regard 
to the use of true quarter or octant elements should be explored numerically. 
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SECTION 5 • 

DESCRIPTION OF THE COMPUTER CODE 


5.1 INTRODUCTION 

Computer code has been prepared to perform all calculations leading to the stiffness 
matrix of Equations(l 14) and (1 1 1) and the stress equations of Equation (116). This 
chapter describes certain features of the code without going into details of the 
programming. 

5.2 ELEMENT GEOMETRY 

The number of nodes on both the exterior and interior boundaries of the hole element 
can be specified by the program user. Because all integrations are performed in the 
first quadrant only, there are the geometric constraints that exterior nodes must be 
placed at the four comers and at the four side midpoints, and interior nodes must be 
placed at 0 = 0®, 90®, 180®, and 270®; see Figure 17. 

The boundaries of the quadrant are divided into sides 1, 2, and 3, as shown in 
Figure 17. It is assumed that the nodes are equally spaced along each side, so that 
the nodal information is input by specifying the number of segments in each side of 
the quadrant; this is done through an array [NSEGJ . For example, in Figure 17, 

[nsegJ =[2 1 3J. 

It is assumed that the hole is centered in the square. The ratio of square dimension, 
2L, to hole diameter, 2R, can be varied as desired by the user; however, it is 
recommended that -^5 2, with possible. 

5.3 NUMBER OF STRESS COEFFICIENTS 

In Section 1 .6 it is shown that there is a relationship between the number of degrees 
of freedom and the number of stress coefficients. What was shown there for a 
general state of deformation must also be true for each of the separate symmetry 
condition states. 
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Figure 17. Hole Element Geometry 
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First/ note that the rigid body motions consist of one SA state (displacement in the 
x-direction)/ one AS state (displacement in the y-direction)/ and one AA state 
(rotation about the center of the element). There is no rigid body SS state. These 
observations establish the values of R to be used when checking the requirement that 
rankCG] = (F-R). 

1 

In the notation of Equation (23), letl_qj^gJ denote all of the nodal displacements on 
the exterior boundaries; clearly, Lq__J- [Oj will suppress all true rigid body motions 
for all symmetry conditions. Then l_q^ Jwill denote the displacements on the hole 
boundary. Let NSEG3 denote the number of segments on the hole boundary (side 3) . 

Then as is shown in Figures 10 through 13, the number of elements in|_q^J (the number 
of degrees of freedom on side 3), denoted by (F-S) in Equation (23), is a function of 
symmetry condition, as follows: 

for SS : (F-S) = NSEG3 + 1 

for SA and AS : (F-S) = NSEG3 (117) 

For AA : (F-S) = NSEG3-1 

In Equations (24) and (25), the array CG^] is the boundary work array associated with 

work done by radial stresses moving through radial displacements on the hole boundary; 

• • 

there will be an array CG.j] for each symmetry condition. Because rank [G^] must 
equal to (F-S). for each symmetry condition, it follows that each Cg'^] must have at 
least (F-S). independent rows. This means there must be at least (F-S). independent 
states of radial stress on the hole boundary. 

2 

Now, consider Equation (64), which shows the assumed form for R r^. Clearly, 
at P = 1 , the stress states associated with ^ and B 2 ^ ore not independent; this is 
also true for Bg^ and B^^ and for many other sets of coefficients. In general, the radial 
stress states which have the same harmonic function of 0 are dependent when evaluated 
along any circle. Therefore, since there must be at least (F-S). independent radial 
stresses on the hole boundary, it follows that there must be a minimum of (F-S). harmonics 
in the stress state for each symmetry condition. 
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For symmetry conditions SS, SA/ and AS, the first two harmonics supply four coefficients; 
for condition AA, the first harmonic supplies two coefficients. These first harmonics are 
never enough because the number of coefficients is not adequate for the number of degrees 
of freedom. Therefore, additional harmonics, supplied by the summations in Equation (64), 
are always necessary. Let AH denote the minimum number of additional harmonics required; 
this is a function of symmetry condition, as follows: 


for SS, SA, AS: AH = (F-S) - 2 
for AA : AH = (F-S) - 1 

Substitution of Equation (117) into Equation (118) gives 

for SS : AH = NSEG3-1 

for SA,AS,AA : AH = NSEG3-2 


(118) 


(119) 


For specified NSEG3, Equation (119) gives the minimum number of additional 
harmonics necessary to provide proper rank to the arrays CG^]. Fewer harmonics 
will introduce spurious rigid body modes involving hole boundary displacements with 
zero displacement on the external boundaries. 

Suppose there are more than (F-S) independent rows in the CG^] array. This is 
perfectly all right as far as rigid body considerations are concerned. The rank of 
array CG^] will equal (F-S), and Equation (25) will be satisfied only by the trivial 
solution = {0}. However, excess independent rows mean excess harmonics in 
the expression for radial stress; and this leads to reduced accuracy in satisfaction of 
any radial traction boundary conditions which might exist on the hole boundary. 

This can be explained with reference to Section 1.7, entitled "Treatment of Boundary 
Conditions." Although the hole boundary has been treated as S^ boundary with 
specified nodal radial displacements, actual traction boundary conditions can be 
approximately satisfied by proper selection of the hole nodal displacements. 

This means there are (F-S) displacements available to provide proper values for 
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radial stress at the hole boundary. If there are more than (F-S) independent stress states 
at the boundary, then there are not enough displacements to permit best satisfaction of 
the traction conditions. If there are exactly (F-S) independent stress states, then in some 
unspecified and probably non-unique manner it is possible to achieve excellent satis- 
faction of traction boundary conditions. This observation, while vague in the details, 

has been numerically confirmed. In open hole test problems described in a following 

1 

section, the radial stresses on the hole boundary were decreased by three orders of 
-2 -5 

magnitude (10 down to 10 or better) by removing excess harmonics from the assumed 
stress state . 


In summary, if the number of additional harmonics is selected to satisfy Equation (119), 
then spurious rigid body modes on the hole boundary will be eliminated and traction 
boundary conditions will be as well satisfied as is possible. With the number of additional 
harmonics known, the total number of stress coefficients depends upon the summation 
limits I and J in Equation (62). It is only necessary to ensure that rankCG]= (F-R) 
for each symmetry condition; this generally requires NSEG3S 3. 

5.4 AREA INTEGRATION 

The internal complementary work leads to the area integral shown in Equation (92) 
and repeated below: 

Ch^] = f Ca^ ( p , 0 )] C c (e)] Ca^ (p,e)]p dpde 
yQuad 


From the structure of Equations (64), (66), and (68), it can be shown that the arrays 
[A^(p, 0)] can be written as 

D1 

[aJ(p,S)]= 2 (0)] 

m = 1 


( 120 ) 


where IJl - I+J+1 

I,J = summation limits shown in Equation (62). 
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In words/ the arrays CA^(p,0)]are decomposed into the summation of products of 
individual powers of p multiplied by functions of 0. 


Substitute Equation (120) into Equation (92) to get 

U\ U1 


2 2/ CC.(6)]C„Al(e):pdpde (121) 

m = l n = l JQuad 


Let P denote the typical integral in Equation (121). Then, it is clear that P has the 
form 

P = / f(p)g(0)dpd0 

yQuad 


( 122 ) 


I ef \ - m+n-2J-5 

where f(p) - p 


(123) 


The integration limits are established with reference to Figure 18, The quadrant is 
divided into two octants identified by the exterior boundary side number. In each 
octant, the upper limit on p is given by 

- I/R 


where 

Define 


max cos <f> 

L/R = ratio of quadrant dimension to hole radius. 


(124) 


F(0) = 


fp 

I max 

ip=l 


f(p) dp 


(125) 


After substituting Equation (123) into Equation (125) the result is 

m+n-2J-4 


1 


F(0) = ' 


m + n - 2J - 4 


\cos <p / 


-1 


, (m+n-2J-4)f^ 


(126) 


In 


l/R 

cos 0 


, (m+n-2J-4) = 0 
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Figure 18. Nol-ation Used for Volume and 

Boundary IntegraHon in fhe Quadrant 
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Now the quadrant integration can be written as the following sum of two octant 
integrations: 



F(0Pg(apdo!^ 



z 

4 


‘^ 0 : 2 ) da 2 


(127) 


where 


0^ = / in octant 1 

7T 

^ 2 ~ 4 ~^2 ' octant 2 


Alternatively the integral is given by 



F(0) C 9 (#) + 9(5 " 0) W 


(128) 


The area integrals in their original forms are too complicated for exact evaluation, 
so numerical integration is necessary. However, it is possible to first reduce the area 
integrals to line integrals through performing the radial integration in the closed form 
of Equation (126). Although numerical techniques are still necessary to evaluate the 
line integrals in Equation (127), the results should be more accurate than what would 
be achieved with comparable effort on the original area integrals. 

5.5 BOUNDARY INTEGRATION 

The external complementary work is calculated from the product of tractions and 
displacements on all three sides of the quadrant. The result is given in Equation (93), 
repeated below. 

The function f(0), used to define the differential boundary length as shown in 
Equation (83), is determined as follows, with reference to Figure 18. 
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Along side 1, the arc length s. is given by s = L tan 0, From this, it follows 

that ds^ = R ^ sec 0 d6; and by comparison with Equation (83), f(6) = (L/R) sec 9. 

If p is introduced from Equation (124), the result can be written as follows: 
max 

2 

f(0)=-^ (129) 

1 

I O 

Alongside 2, $2 = LCl-tan(^ - 0)] . Therefore, ds 2 = R'^ sec (•^ - 0)d0. Since 

l/r 

p = , it follows once again that f(0) is given by Equation (129). 

max /IT 

cos(^ - 0) 

Along side 3, Sg = R0 and ds^ = Rd0 . 

Therefore, f(0) = 1 . 

The tractions on sides 1 and 2 are in the x and y directions while the stresses are 
assumed in polar coordinates. Therefore, at each point on the exterior boundary, 
the stresses are first evaluated and then transformed to cartesian components by well- 
known equations such as given in [l^ p.67] . 

On side 3, the radial traction Is immediately given by the negative of the radial stress 
component. 

The integrals in Equation (93) are finally evaluated numerically as discussed in the 
next section. 

5.6 NUMERICAL INTEGRATION 

In the numerical evaluation of the area integrals, each octant, described by 45°, 

is divided into a user specified number of integration intervals. That is, each octant is 
divided into a specified number of equal Increments of the angle variable. Then the 
integral in Equation (127) is evaluated in each interval by Gaussian quadrature with a 
user specified number of Gauss points. 
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For the boundary integration, each side is first divided into segments by the nodes. Then 
each segment is further subdivided into integration intervals, with the number of intervals 
possibly different for each side. Finally, the integral is evaluated in each interval by 
Gaussian quadrature. 

As mentioned above, the number of Gaussian quadrature points per interval is user 
specified. However, the existing computer code limits the choices to 3, 5, 7, or 9 
points. This choice can, of course, be expanded if desired. 


76 


SECTION 6 
SAMPLE PROBLEMS 


6J INTRODUCTION 

NASTRAN was chosen as the host program for this element because it is such a widely used, 
well-known program with a large element library and many capabilities. The hole element 
implemented into NASTRAN is a variab le-noded total element and is defined via ADUM7, 
CDUM7, and PDUM7 input cards. The format of these cards are given in Appendix A and 
the default geometry of the element is shown in Figure 19. The user may specify the 
number of nodes on each of the three boundaries shown in the figure. The nodes on the 
exterior of the element have two cartesian degrees of freedom whereas the nodes on the 
hole wall have only a single radial degree of freedom. Tangential degrees of freedom at 
the hole wall do not appear in the element fomiulation because the stiffness matrix was 
derived through a complementary virtual work formulation based on assumed stresses which 
satisfy the tangential traction-free conditions on the hole boundary (see Section 2.1). The 
ratio of element diameter to hole diameter is determined from grid point locations and is 
thus a user choice; however, a ratio of four is recommended. The material properties may 
be isotropic or orthotropic and are defined by MATl and MAT2 cards. Characteristics 
which are not defined for this element include mass properties, thermal loading, material 
non-linearity, and differential stiffness. 

The example problems which follow were solved using the NASTRAN program at Lockheed- 
Georgia Company with the orthotropic hole element implemented. They are intended to 
demonstrate the performance of the element In a variety of applications. 

6.2 ISOTROPIC PLATE WITH UNIAXIAL TENSION 

The orthotropic hole element can be used in isotropic applications by simply specifying 
isotropic elastic coefficients. This is demonstrated by presenting the solution for an infinite 
isotropic plate with unit uniaxial tension at infinity (see Figure 20). The finite element 
model is shown in Figure 21 . The plate width is 20 times the hole diameter, which repre- 
sents an essentially infinite plate. Poisson's ratio is set at 0.3205. 
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Figure 20. Unif Tensile TracHon in the Y-Direction 
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Table I presents results for the three stress components at four values of radius and four 
values of 0. The hole element stiffness and stress matrices were based on array 
[nsegJ = [l 1 sj with radial variation given by powers of p from p ^ to p^^ . 


The hole boundary is given by r= 1 . Note the excellent satisfaction of the zero radial 
stress boundary condition; the results were typically 10 ^ times the magnitude of external 
loading. The circumferential stress results are also very good, with an error of only 0.2% 
on the maximum tension stress concentration factor and 1 .5% on the maximum compression 
stress. The zero shear stress condition is of course perfectly satisfied since it was included 
in the original stress assumptions. 


The stresses away from the hole boundary are also very well represented by the hole 
element results. 


6.3 ISOTROPIC PLATE WITH HALF COSINE LOADED HOLE 

A load transfer problem for which a comparison solution is available [3] is that of an 
infinite isotropic plate with Poisson's ratio equal to 0.25, loaded with uniform traction 
on one end and a cosine variation of radial tractions on the hole boundary (see Figure 22). 
The total load on the end is equal to 2P, and the radial tractions on the loaded hole 
boundary are given by 


4P 

T = — COS0, 
r ff 



The finite element model is shown in Figure 21; and Table II presents results in the region 
around the hole, based on [_NSEgJ “ 1 and radial terms from p to p . 

The value of P is set equal to 1 .0. 


The stresses at the hole boundary are very well represented by the finite element results. 

The radial stresses follow the cosine curve over the loaded region and are of the order of 
—9 

lO”"^ on the unloaded portion of the boundary. The variation of circumferential stress 
is reasonably accurate, with an error of 3% on the maximum stress concentration. 

In the region around the hole, the results from the hole element solution continue to 
give a good description of the stress state. 
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TABLE II 

STRESSES FOR AN ISOTROPIC PLATE WITH LOAD TRANSFER 


r = 1.2 


r = 1 .5 


.271 
.099 
-0.649 
-0.041 
- 0.012 
0.003 
0.003 


-0.995 

-0.844 

-0.445 

0.038 

0.055 

0.023 

0.006 


-0.726 

-0.604 

-0.288 

0.045 

0.094 

0.055 

0.035 


-1.273 

-1.103 

-0.637 

0 

0 

0 

0 


-0.977 

-0,829 

-0.429 

0.034 

0.078 

0.038 

0.024 


-0.726 

-0.609 

-0.293 

0.037 

0.086 

0.049 

0.033 


’*0 

* ★ 

T * 

^r0 


0.372 

0.412 

0 

0 

0.467 

0.470 

0 

0 

0.659 

0.612 

0 

0 

0.834 

0.810 

0 

0 

0.415 

0.373 

0 

0 

0.056 

0.053 

0 

0 

-0.094 

-0.065 

0 

0 


0.348 

0 

0 

■iii 

0.381 

0.035 

0.021 

0.460 

0.443 

0.046 

0.029 

0.477 

0.415 

-0.032 

-0.031 

0.289 

0.274 

-0.097 

-0.081 

0.070 

0.088 

-0.064 

-0.051 

-0.015 

-0.011 

0 

0 


0.241 

0 

0 

mm 

0.253 

0.042 

0.026 

0.295 

0.266 

0.046 

0.028 

0.266 

0.231 

-0.047 

-0.045 

0.179 

0.152 

-0.119 

-0.105 

0.062 

0.055 

-0.083 

-0.070 

0.010 

0.012 

• 0 

0 


Finite Element Results 
Reference 3 





In the finite element model of the structure, the tractions ore applied as concentrated 
nodal forces. On the hole boundary, these nodal forces are the work equivalent loads 
based on the cosine variation of the radial traction between nodes and the assumed 
radial displacement of Equation (79). The results are as follows. 

Let A = angular increment subtended by arc length between two nodes 

on the hole boundary 



These results can be partially checked as follows. Consider a hole boundary with N 
segments per quadrant. Then the summation of nodal force components in the y 
direction is 

y' F = ~ A + 2 (— A cos^ a)+ 2(— A cos^2a) +. . .+2 [— A cos^(N-l) A] 
y IT It TT IT 

= — aCI +2(cos^ A+ cos^ 2a +. . .+ cos^ (N-1) A ] 

y TT 

Since A = , it follows that 
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It can be shown that 


N-1 

y 2. 

L-t cos I 


i = 1 



N-1 
2 ' 


N 52 


Therefore, ^ [l + 2 (^)]= 2P 

as required for equilibrium. 


6.4 ±45° LAMINATED PLATE 

The example problems with isotropic material serve to confirm the basic theoretical 
approach and to partially validate the computer solution. There remains only to 
demonstrate satisfactory performance in the intended applications with orthotropic 
materials. The first material selected has material properties as follows; 

E =E =E, G = 1 .697528E, v =v =0.735 
X y xy yx xy 

Note that even though E^ = E^, this material is rather "far" from isotropic, with some 
unspecified measure, due to the large value of Poisson's ratio and the fact that E,G, and i> 
do not satisfy the isotropic relationship. These properties describe the ±45“ graphite/ 
epoxy laminate considered in [4]. 

The first loading is unit tensile tractions in the y-direction, applied as shown in Figure 20 
with the finite element model given in Figure 21 . The traction boundary conditions at 
the hole boundary are satisfied, with radial stress of the order 10 ^ times the external 
loading. Figure 23 shows the variation of circumferential stress around the hole boundary 
for [nSEgJ = jj 1 ^ and two different sets of powers of p. Also shown is the exact 

solution for an infinite plate, given by 

_ 3 . 05786 cos^ 0 - 1 

® 4 2 

2 . 8809072(cos 0 - cos 0) + 1 
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which can be deduced from [5, p.l75]. The orthotropic hole element results contain 
the essential feature that the maximum stress concentration does not occur at the 
familiar location of 0 = 0®, The worst solution shown has an error in maximum stress 
of only 3.8%, and the best error is only 0.7%, 

The hole loaded with half-cosine distributed tractions has also been solved; see Figure 22. 
The resulting radial stresses on the boundary are essentially identical to the results achieved 
in the case of isotropic material, as should be. The maximum circumferential stress on 
the hole boundary is Tg = 0.79 at 0= 125®, which can be compared with the hole element 
results for isotropic material of T« = 0.84 at 0= 85®. Again, it is seen that the location 
of maximum stress has shifted around the circle from the familiar location. 


6.5 PLYWOOD 

The next material selected has properties as follows (see Section 2.5 for notation): 


E =E, E =2E, G =0.11667E, v =0.036 
y X xy xy 


These properties describe a plywood material for which [5] contains many results for an 
infinite plate. Note that the plywood properties are quite different from the ±45® 
lamination properties, both in magnitudes and relationships. This provides further checks 
on the ability of the hole element to handle different types of orthotropic material. 


The types of loading considered are unit tensile tractions in the y-direction (Figure 20), 
unit compressive tractions in the x-direction with constraint in the y-direction (Figure 24a), 
and uniform unit internal pressure on the hole boundary (Figure 24b). For all cases, the 
finite element model is shown in Figure 21; and [nSEC^=[j 1 dj , with powers of 


r -5 , +1 

p from p to p , 


For the uniaxial tension, there is the usual excellent satisfaction of traction boundary 

-6 


condition, with T of the order 10 
r 


Circumferential stresses around the hole boundary 
are shown in Figure 25. Also shown are the exact results for an infinite plate, given by 
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b) Interna I Pressure 


Figure 24, External Loadings 
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Figure 25. Circumferential Stress at the Open Hole Boundary for Plywood 
with Unit Uniaxial Tensile Traction 


which can be deduced from [5, p.l75]. The general features of the variation of stress 
are demonstrated by the hole element results. There is an 8% error in maximum stress. 


The constrained compression results satisfy the hole boundary conditions, with T of 

-5 

the order 10 . Circumferential stresses around the boundary are shown in Figure 26. 

The exact results for infinite plate are given by 


^0 = 


1.1155 - 6.5174 sin 0 


-13.9989 sin'^0 + 12.9989 sin^0 + 2 


deduced. from [5, p.l79]. There is an 18% error in maximum stress. 


The internal pressure results provided boundary radial stresses of the magnitude -1 .00002 
which is again excellent satisfaction of boundary conditions. Circumferential stresses 
around the boundary are shown in Figure 27, compared to infinite plate exact results 
given by [5, p.l73] 


r 


0 


4 


8829 - 15.8432 sin^0 + 13.9989 sin'^0 
2 + 12.9989 sin^0- 13.9989 sin'^0 


The general shape of the distribution is indicated, but there is a 21% error in the 
maximum stress. Also, there is a small region where the stress becomes compression, 
which is not to be expected for the internal pressure loading and is not shown by the 
exact solution. 


6.6 MULTIPLE HOLES 


Figure 28 shows a finite element model of a plate containing five holes. For economy 
in modeling, it is desirable to develop the hole element stiffness matrix only once and 




0 30 60 90 

$ - deg. 


Figure 26. Circumferential Stress at the Open Hole Boundary for Plywood 
with Unit Constrained Compression 



Figure 27. Circumferential Stress at the Hole Boundary for Plywood 
with Unit Internal Pressure on the Boundary 



Figure 28. Finite Element Model for o Plate with 
Multiple Holes 
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then transform this matrix as necessary to represent all the hole elements. This is possible 
in NASTRAN by using the congruency feature, and this feature has been checked by 
determining the solution to Figure 28 with tension applied parallel to the row of holes. 
The material is isotropic with Poisson's ratio of 0.3205. 


The maximum stress concentration factor is given in [6, p.94] for an intermediate hole in 

an infinite row of holes. In order to simulate the infinite row, the displacements at each 

loaded end of the model were constrained to be uniform across the width of the plate. 

The finite element result for stress concentration is 2.78, based on [_NSEC^ = |j 1 sj 

-5+1 

and powers of p from p to p . Reference 6 gives 2,65; there is a 5% difference in 
results. 


6.7 EFFECT OF VARYING THE HOLE SIZE 

All the preceding examples used a hole element with the ratio of square dimension to 
hole diameter equal to four, which is the recommended value. However, the computer 
code does allow the user to vary this ratio, and this feature of the code was tested. 

A problem for which theoretical results are available is that of an infinite strip under 

unit uniaxial tensile traction with a hole diameter equal to half the plate width. The 

finite element model consists of a hole element which spans the strip width plus sufficient 

linear strain quadrilaterals added symmetrically to achieve a model with length to width 

ratio of five. The hole element stiffness matrix is based on (_NSEgJ = (_2 2 4j , 

-5 +1 

with powers of p from p to p . The material is isotropic, with Poisson's ratio of 
0.25. 

The maximum stress occurs on the hole boundary at 6 = 0°. Reference 5, page 84, 
gives a stress concentration factor of 4.32, The finite element result is 3.95, for an 
8% difference; this is an excellent result for such a relatively crude model. 
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SECTION 7 

IMPLEMENTING THE HOLE ELEMENT INTO NASTRAN 


The orthotropic hole element has been implemented into NASTRAN through the dummy 
element capability. This capability permits the user to enter his own element sub- 
routines for the purpose of generating the stiffness matrix contributions and for the 
computation and output of stresses [7, section 6.8.5]. 

The hole element has been implemented as a DUM7 element. The procedure used to 
implement the element is as follows. 

• Create a subroutine KDUM7 which computes and outputs to functional 
module EMG the element stiffness matrix. Subroutine KDUM7 calls 15 
new auxiliary subroutines in performing its task. Table III gives the 
function of each of these subroutines. Subroutine EMGPRO was updated 
to increase the allowable number of words in the element summary table. 
LINKS was remapped and the overlay structure was modified. 

• Create two subroutines SDUM71 and SDUM72 to compute and output to 
functional module SDR2 the element stresses. Subroutine SDUM72 calls 
one new auxiliary subroutine, STRSIS, which evaluates the elements of 
the assumed stress state. Subroutine SDR2E was changed to call SDUM72 
once for each line of stress output. This was required because this element 
produces many lines of stress output for each element in the model. LINK13 
was remapped. 

• Subroutine GPTABD was modified to increase the number of words SDR2 
passes from phase 1 element stress recovery routines to phase 2 routines. 

The following subroutines were changed to increase the allowable number 
of data items on the CDUM7 bulk data card: 
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TABLE III 

NEW SUBROUTINES WRITTEN FOR NASTRAN 


NAME 

FUNCTION 

ORTHO L 

Computes the orthottopic hole element stiffness matrix 

SETRUN 

Sets up dimension blocks and controls the execution 

STRESl 

Gives stresses for shear free hole boundary 

BNDWRK 

Performs exterior boundary integration around a quadrant 

BNDTRC 

Evaluates tractions on outer boundary 

RDSTRS 

Evaluates radial stress on inner hole boundary 

VOLINT 

Evaluates the volume integrals in a quadrant 

INTGND 

Forms the volume integrand 

ELAS 

Evaluates elastic property matrices 

VOLSTl 

Evaluates stress terms in the volume for shear free hole 

INVS 

Inverts symmetric square array 

AMTRXl 

Evaluates nodal transformation matrices for total element 

EQUIL 

Checks total element equilibrium 

HLCOND 

Modifies arrays to account for open holes 

SQCHK 

Evaluates stresses in quadrant of open hole element 


IFX4BD 

IFPDCO 

IFS2P 

IFS5P 

IFX7BD 

IFP3B 

IFS3P 

IFXDBD 

IFP 

IFSIP 

IFS4P 

IFX5BD 


LINKl was remapped . 

• The following subroutines were modified to account for the increased 
size of the element connectivity table, ECT; 

GPl TAIA TAIETD 

GP2 TAIB TAIH 

LINK2 was remapped. 

Listings of all new and modified subroutines are supplied under separate cover. 
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APPENDIX A 

NASTRAN INPUT CARDS FOR THE ORTHOTROPIC 
HOLE ELEMENT 


A-1 



BULK DATA DECK 


Input Data Card ADUM7 Hole Element 


Description; Defines attributes for the N-node hole element 


Format and Example; 



NG 

NC 

NP 

ND 

fS^ 






24 

29 

11 

3 







Field; Contents 


NG Number of grid points connected by CDUM7 element 
(Integer S 64) 

NC Number of additional entries on the CDUM7 connection card. All entries which 

follow the grid point identification numbers are additional. (Integer 5 0) 

NP Number of additional entries on the PDUM7 property card. All entries which follow 

the MID entry are additional. (Integer 5 1) 


ND Number of displacement components at each grid point used in generation of 

differential stiffness matrix (Always input 3) 

Remarks; I . The value of NC given in the example above is not typical but is consistent 
with the CDUM7 and PDUM7 examples. Typically, if additional entries are 
present on the CDUM7 card, they will consist of only the NS EG array (i.e. NC=3). 
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BULK DATA DECK 


Input Data Card CDUM7 Hole Element Connection 

Description; Defines an N-node hole element of the structural model 


Format and Example: 


CDUM7 

EID 

PID 

G1 

G2 

G3 

G4 

-etc- 

GN 

abc 

CDUM7 

101 

10 

12 

15 

96 

97 


11 

ABC 


4bc 

NSEGl 

NSEG2 

NSEG3 

NINTl 

NINT2 

NINT3 

NIPB 

NIO 

def 

+BC 

1 

1 

4 

2 

2 

2 

5 

2 

DEF 


+ef 

NIPV 

LPl 

LP2 

LP3 

LP4 

LNl 

LN2 

LN3 

ghi 

4€F 

5 

3 

3 

3 

3 

3 

3 

3 

GHI 


+hi 

LN4 

ITOTRT 

IQDTRT 

ICOND 

lACHK 

lELMNT 

IBNDRT 

IVOLRT 

iki 

4HI 

3 

0 

0 

0 

0 

0 

0 

0 

JKL 


mm 

IHINVS 

ISCl 

ISC2 

1SC3 

NLC 





+KL 

0 

1 

4 

1 

4 






Field 

EID 

PID 

G1....GN 


Contents 

Element identification number (Integer > 0) 

Identification number of a PDUM7 property card (Integer > 0) 

Grid point identification numbers of connection points 
(Integer > 0, G1 ^ G2. . . GN) 
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BULK DATA DECK 


Input Data Card CDUM7 Hole Element Connection (Continued) 


Field 

NSEGi 

NINTi 

NIPB 

NIO 

NIPV 

LPi 

LNi 

ITOTRT 

IQDTRT 

ICOND 

lACHK 

IEUV\NT 

IBNDRT 


Contents 

Number of segments on boundary i (see Figure A-1, Integers, Defaults 
are NSEGI = 1, NSEG2 = 1, NSEG3 = 4) 

Number of integration intervals per segment for boundary i 
(Integer, Defaults = 2) 

Number of Gaussian integration points per interval for boundary 
integration (Integer, Default = 5) 

Number of integration intervals per volume octant 
(Integer, Default = 2) 

Number of Gaussian integration points for volume integral 
(Integer, Default = 5) 

Number of stress terms with powers of r > - 2. Input for each symmetry 
condition where i = 1 is symmetric-symmetric, i = 2 is symmetric- 
antisymmetric, i = 3 is antisymmetric-symmetric, and i = 4 is antisymmetric- 
antisymmetric. (Integer^ 3, Defaults = 3) 

Number of stress terms with powers of r < -2. Input rules are the same as 
for LPi. (Integers 3, Defaults = 3) 

If greater than zero, causes a print of the elemental stiffness and stress 
matrices. (Integer, Default = 0) 

If greater than zero, causes a print of intermediate stiffness and stress 
matrices for a quadrant of the element (Integer, Default = 0) 

If greater than zero, causes a static condensation of the element stiffness 
and stress matrices through elimination of the degrees of freedom on a 
traction free hole wall (Integer, Default = 0) 

If greater than zero, causes a print of the transformation matrix which 
relates total and quadrant displacements (Integer, Default = 0) 

When equal to 0, produces total element stiffness and stress matrices. 

If -1, quadrant element matrices are generated, and if +1 , half element 
matrices are generated. (The only option currently available is 0 ) 

If greater than zero, prints the boundary work matrix (Integer, Default = 0) 

A-4 



BULK DATA DECK 


Input Data Card CDUM7 Hole Element Connection (Concluded) 

Field Contents 

IVOLRT If greater than zero, prints the volume work matrix (Integer, Default = 0) 

IHINVS If greater than zero, prints the identity matrix resulting from the check 

of the inversion of the volume work matrix (Integer, Default = 0) 

If lELMNT is +1 or -1, symmetry conditions 
must be defined. In the program code, symmetry 
conditions are set with DO loops such as 
DO XX 1 = ISCl, ISC2, ISC3 

See the description of LPi for definition of symmetry conditions. 

The user should chose ISCl, ISC2, ISC3 to implement the conditions 
he wishes. NLC is the number of times the DO loop is executed. 

This is for dimension purposes. (Integer, Defaults; ISCl = 1, ISC2 = 4, 
ISC3 = 1, NLC =4) 

Remarks: 1 . Element identification numbers must be unique with respect to £ll 

other element identification numbers. 

2. The element must be planar with a square outer boundary and a 
circular inner boundary. The nodes on the outer boundary, as well 
as the inner boundary, must be equally spaced. 

3. Nodes on the hole wall have only radial degrees of freedom, therefore 
tangential restraints are required. This is necessary because in a general 
finite element model there are six degrees of freedom associated with each 
node and the stiffness matrix will be singular if restraints are not applied to 
those degrees of freedom not included in the element formulation. The 
tangential restraints are imposed only to eliminate the singularity in the 
stiffness matrix. The actual tangential displacements are not zero, and 
circumferential strain is developed at the hole boundary. 

4. For reference purposes, the values given for all additional entries on the 
example input card are default values. If default values are acceptable 
for all additional entries, they may be omitted from the input. If a 
variable is given a non-default value, all preceding variables must be 
defined; all subsequent variables may be omitted. 


ISCl, 

1SC2, 

ISC3 

NLC 
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BULK DATA DECK 


Input Data Card PDUM7 Hole Element Property 


Description: Defines the properties of a N-node hole element. 
Referenced by the CDUM7 card. 


Format and Example: 


PDUM7 

PID 

MID 

T 

DTH 

THMX 

NR 

RI 

R2 

abc 

PDUM7 

10 

2 

.1 

10. 

359. 

4 

1. 

1.5 

ABC 


•tbc 

R3 

-etc- 

RN 

PRTl 

PRT2 

PRT3 




+BC 

2.0 


3.5 

0 

0 

0 





Field: 


Contents 


PID 

MID 

T 

DTH,THMX, 

NR,Ri 


PRTl 


.Property identification number (Integer > 0) 

Material identification number (Integer > 0) 

Thickness (Real) 

Defines stress output locations. Stresses will be computed and printed 
for the NR radii defined by R1 ,R2,,..,RN and at angles beginning at 
0° and incremented by DTH, up to a maximum angle of THM.X 
(See Figure A-2). DTH, THMX, and Ri are real, NR is integer S20. 
Defaults: DTH = 10®, THMX = 359.999, NR = 4, Ri equally spaced 
from hole wall to element boundary. Note: Ri are normalized to a 
unit radius. 

If greater than zero, causes a print of the total element stress matrices 
( Integer) 


PRT2 

PRT3 


If greater than zero, causes a print of displacements in the element 
local coordinate system (Integer) 

If greater than zero, causes a print of the stress coefficients associated 
with elemental nodal displacements (Integer) 
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Remarks: 


1 . 


All PDUM7 cards must have unique identification numbers. 


2. Stress output for this element (see Figure A-3) may be interpreted 
as follows: 

Field SI - The radius, r, at which stresses are computed. 

Field S2 - The angle, 0, at which the stresses in fields 
S3, S4, and S5 are computed. 

Field S3 - Tj, 

Field S4 - Tg 

Field S5 - 

Field S6 - The angle, & at which the stresses in fields 
S7, S8, and S9 are computed 

Field S7 - 
Field S8 - 
Field S9 - 

Note that each line of output has stress values for two 
values of 6. 
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Figure A-3. Sample NASTRAN Stress Output for the CDUM7 Element 



APPENDIX B 

INTERPRETATION OF OUTPUT 


Figure B-1 gives on example of the output which permits the user to confirm that the 
hole element input was correct. Number of segments per side shows the entries in array 
(_NSEGJ / and number of intervals per segment shows the entries in array |_NINTj . 
Number of terms with powers of R > -2 (R < -2) gives the values of II , 12, 13, and 14 
(Jl, J2, J3, and J4) in Equation (63). For each symmetry condition, number of stress 
coefficients gives the correct number based on the minimum number of harmonics necessary 
to eliminate spurious rigid body motions; these numbers are determined by the program 
based on input of NSEG3 and values of I and J in Equation (63). 


The material properties are described first by the values of E , E , y , and G . The 
number of elastic constants, NEC, displays the number of independent constants; NEC = 2 
denotes an isotropic material. Constants Cl , C2, C3, and C4 are the coefficients defined 
by Equation (39); and CBAR denotes the C defined in Equation (51). 


The element geometry is described by the value of. the hole radius, the L/R value 
(see Figure 17), and the plate thickness. 


The quantities ISCl, ISC2, and ISC3 are used by the program to determine the symmetry 
conditions involved in the derivation of the hole element stiffness matrix; NLC denotes 
the number of load conditions used. In the present form of the program, which incorporates 
only the total element, ISCl, ISC2, ISC3, and NLC must have their default values 
of 1,4, 1 , and 4 respectively. 

NBMAX is simply the largest element in the array of stress coefficients. MXLPNl denotes 
the largest number of powers of p for all symmetry conditions; this is the maximum IJl in 
Equation (120). The values NBMAX and MXLPNl are used for DIMENSION purposes in 
the program. 

The displacement degrees of freedom are denoted as follows: 
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Figure B-1 . Program Output Concerning the Hole Element 





EOUIL 
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Figure B-2. Equilibrium Check for Total Element 
Stiffness Matrix 
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NQDOF 


= number of degrees of freedom for fhe quadrant element 

NQHDOF = number of degrees of freedom on the hole boundary of 
the quadrant 

NQEDOF = number of degrees of freedom on the outer boundaries 
of the quadrant 

NTDOF = number of degrees of freedom for the total element 

When executing a NASTRAN solution, there is a certain amount of computer core 
available for the orthotropic hole element code. If the available space is sufficient, 
the value of JMAX will give the number of memory locations necessary to provide for 
all the arrays used in the hole program. If the available space is not sufficient, a 
message will appear which gives the minimum number of additional storage locations 
which must be made available through the NASTRAN HICORE specification (it is 
recommended that this minimum number be increased by at least 2000 to supply a margin 
of safety); in this case, program execution is terminated. 

The remaining information shown in Figure B-1 consists of the number of Gaussian 
integration points per interval for the boundary integration (NIPB), number of integration 
intervals per volume octant (NIO), and number of integration points per interval for the 
volume integration (NIPV). 

Figure B-2 gives a typical final equilibrium check for the stiffness matrix of the total 
orthotropic hole element. The values have been nonnalized with respect to the diagonal 
element in each column of the stiffness matrix. 

Figures B-1 and B-2 give the minimum hole element output currently furnished to the user. 
However, much other output can be obtained if desired by exercising the WRITE options 
discussed in the description of the CDUM7 data card. These options are further described 
as follows: 
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ITOTRT > 0 


IQDTRT > 0 


will print the total stiffness matrix [K] in Equation (111) 
and the total stress matrices RCS^], with CS^] given in 
Equation (115); 

will print the quadrant stiffness matrices [KQ^] in 
Equation (100) and the quadrant stress matrices RCBQP ], 
with CBQ^] given in Equation (94); 


ICOND > 0 


lACHK > 0 
IBNDRT>0 
IVOLRT > 0 
IHINVS > 0 


will result in transformed stiffness and stress matrices for 
quadrant elements with traction free hole boundaries; output 
will consist of intermediate steps in the transformation, the 
final modified stress and stiffness matrices, and stresses on 
the hole boundary due to unit values of the nodal displace- 
ments on the outer boundaries of the element, 

will print the arrays CA^ 1 in Equation (109); 

will print the boundary work arrays [G^ ] in Equation (93); 

will print the volume work arrays CH^ ] in Equation (92); 

i * “1 

will print the products CH^] CH^] and provide some check 
on the accuracy of the inverse of CH^]. 


This completes the description of the output which is associated with the development of 
the hole element stiffness and stress matrices. Final output, consisting of stress coefficients 
and values of stress components at specified locations, is discussed in the description of 
the PDUM7 data card. 
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